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^ ■ Abstract 
■ 

^ We investigate the coherent propagation of dilute atomic Bose-Einstein condensates 

^ . through irregularly shaped billiard geometries that are attached to uniform incoming and 
^ i outgoing waveguides. Using the mean-field description based on the nonlinear Gross-Pita- 
Q-f evskii equation, we develop a diagrammatic theory for the self-consistent stationary scatter- 
ing state of the interacting condensate, which is combined with the semiclassical represen- 
tation of the single-particle Green function in terms of chaotic classical trajectories within 
the billiard. This analytical approach predicts a universal dephasing of weak localization 
in the presence of a small interaction strength between the atoms, which is found to be 
in good agreement with the numerically computed reflection and transmission probabilities 
of the propagating condensate. The numerical simulation of this quasi- stationary scatter- 
ing process indicates that this interaction-induced dephasing mechanism may give rise to 
a signature of weak antilocalization, which we attribute to the influence of non-universal 
short-path contributions. 
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1. Introduction 

Recent technological advances in the manipulation of ultracold atoms on microscopic 
length scales have paved the way toward the exploration of scattering and transport phe- 
nomena with coherent interacting matter waves. Key experiments in this context include 
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the creation of flexible waveguide geometries with optical dipole beams [H and on atom 
chips jl, [sj , the coherent propagation of Bose- Einstein condensed atoms in such waveguides 
by means of guided atom lasers [3, H, H, 0| , the realization of optical billiard confinements 
[sl, 0, US] and microscopic scattering and disorder potentials for cold atoms [ll|, 12|, as well 
as the detection of individual atoms within a condensate through photoionization on an 



atom chip [13|. Moreover, it was recently demonstrated [ij] that artificial gauge potentials 
can be induced for cold atoms, which lead to a breaking of time-reversal invariance in the 
same way as do magnetic fields for electrons. Such artificial gauge potentials can, e.g., be 
implemented by means of Raman dressing with two laser beams that include a finite orbital 
angular momentum [l^, 16, 17]. Together with the possibility of combining different atomic 
(bosonic and fermionic) species and of manipulating their interaction through Feshbach res- 
onances, the combination of these tools gives rise to a number of possible scattering and 
transport scenarios that are now ready for experimental investigation. 

A particularly prominent quantum transport phenomenon in mesoscopic physics is weak 



localization [18|, ll9| . This concept refers to an appreciable enhancement of the reflection (or. 



in the solid-state context, of the electronic resistance) in the presence of a two- or three- 
dimensional ballistic or disordered scattering region, as compared to the expectation based 
on a classical, i.e. incoherent, transport process. This enhancement, which in turn implies a 
reduction of the transmission (or of the electronic conductance) due to current conservation, 
is in particular caused by "coherent backscattering" , i.e. by the constructive interference 
between backscattered classical paths and their time-reversed counterparts, which was first 
observed in experiments on the scattering of laser light from disordered media 20|, |2l|]. In 



the solid-state context, weak localization is most conveniently detected by measuring the 
electronic conductance in dependence of a weak magnetic field that is oriented perpendicular 
to the scattering region, such that it causes a dephasing between backscattered paths and 
their time-reversed counterparts. A characteristic peak structure at zero magnetic field is 



then typically observed [22|, |23 



From the electronic point of view, the presence of interaction between the particles that 
participate at this scattering process is generally expected to give rise to an additional 



dephasing mechanism of this subtle interference phenomenon [2J, |25|, l26|. In the context 
of ultracold bosonic atoms, this expectation is partly confirmed by previous theoretical 
studies on the coherent propagation of an interacting Bose-Einstein condensate through 
a two-dimensional disorder potential 
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which employed numerical simulations as well 
as diagrammatic representations based on the mean-field description of the condensate in 
terms of the nonlinear Gross-Pitaevskii equation. This study did indeed reveal a reduction of 
the height of the coherent backscattering peak with increasing effective interaction strength 
between the atoms. It also predicted, however, that this coherent backscattering peak might 



turn into a dip at finite (but still rather small) interaction strengths 27|]. This scenario is 



reminiscent of weak antilocalization due to spin-orbit interaction, which was observed in 
mesoscopic magnetotransport [28]. 

In order to gain a new perspective on this novel phenomenon, we investigate, in this 
work, the coherent propagation of Bose-Einstein condensates through ballistic scattering 
geometries that exhibit chaotic classical dynamics. Such propagation processes can be ex- 
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perimentally realized by guided atom lasers in which the optical waveguides are locally "de- 
formed" by means of additional optical potentials, e.g. by focusing a red-detuned laser from 
a different direction onto this waveguide as was done in the experiment of Ref. 0]. Alterna- 
tively, atom chips j2l, Isj or atom-optical billiards js], [q], [loj could be used in order to engineer 
chaotic scattering geometries for ultracold atoms. From the theoretical point of view, the 
wave transport through such scattering geometries can be described using the semiclassical 
representation of the Green function in terms of classical trajectories. The constructive in- 
terference of reflected trajectories with their time-reversed counterparts gives then rise to 



coherent backscattering [29|, while a complete understanding of weak localization, in partic- 



ular the corresponding reduction of the transmitted current, requires additional, classically 



correlated trajectory pairs [30|, |31 



In order to account for the presence of atom-atom interaction on the mean-field level 
of the nonlinear Gross-Pitaevskii equation, we combine, in this paper, the semiclassical ap- 
proach with the framework of nonlinear diagrammatic theory developed in Refs. 32|, |33|, l34 



For the sake of simplicity, we shall, as is described in Section |21 restrict ourselves to ideal 
chaotic billiard dynamics consisting of free motion that is confined by hard-wall boundaries. 
Since such billiard geometries give rise to uniform average densities within the scattering 
region, we can, as demonstrated in Sections [3] and HI derive explicit analytical expressions for 
the retro-reflection and transmission probabilities as a function of the effective interaction 
strength. As shown in Section these expressions agree very well with the numerically 
computed retro-refiection and transmission probabilities for two exemplary billiard geome- 
tries as far as the deviation from the case of noninter acting (single-particle) transport is 
concerned. On the absolute scale, however, the height of the weak localization peak is re- 
duced in this noninteracting case by the presence of short-path contributions, in particular 
by self-retracing trajectories, which, as shown in Section |5l consequently turn this peak into 
a finite dip in the presence of a small interaction strength. We shall therefore argue in Sec- 
tion [H] that such short-path contributions are at the origin of this weak antilocalization-like 
phenomenon. 



2. Setup of the nonlinear scattering process 

We consider the quasi-stationary transport of coherent bosonic matter waves through 
two-dimensional waveguide structures that are perturbed by the presence of a wide quantum- 
dot-like scattering potential. Such propagating matter waves can be generated by means of 
a guided atom laser 0, IH where ultracold atoms are coherently outcoupled from a trapping 
potential that contains a Bose-Einstein condensate. The control of the outcoupling process, 
which, e.g., can be achieved by applying a radiofrequency field that flips the spin of the 
atoms in the (magnetic) trap permits one, in principle, to generate an energetically 
well-defined beam of atoms that propagate along the (horizontally oriented) waveguide in 
its transverse ground mode [iij . This waveguide, as well as the quantum-dot-like scattering 
potential, can be engineered by means of focused red-detuned laser beams which provide 
an attractive effective potential for the atoms that is proportional to their intensity. The 
restriction to two spatial dimensions can, furthermore, be realized by applying, in addition. 
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a tight one- dimensional optical lattice perpendicular to the waveguide (i.e. oriented along 
the vertical direction). 

The central object of study in this work is the phenomenon of weak localization. In 
the context of electronic mesoscopic physics, this quantum interference phenomenon can be 
detected by measuring the electronic conductance, which is directly related to the quantum 
transmission through the Landauer-Biittiker theory 36|, |37|, |38|] , as a function of the strength 
of an externally applied magnetic field which breaks time-reversal invariance within the scat- 
tering region. Such a time-reversal breaking mechanism can also be induced for cold atoms 



14 Il5l . Il6l . llTj . e.g., by coherently coupling two intra-atomic levels via a STIRAP process 



using two laser beams of which one involves a nonvanishing orbital angular momentum [15 



This gives rise to an effective vector potential in the kinetic term of the Schrodinger equation, 
which is assumed such that it generates an effective "magnetic field" that is homogeneous 
within the scattering region and vanishes within the attached waveguides. 

The main purpose of this study is to investigate how the scenario of weak localization 
is affected by the presence of a weak atom-atom interaction within the matter-wave beam. 
In lowest order in the interaction strength, the presence of such an atom-atom interaction 
is accounted for by a nonlinear contribution to the effective potential in the Schrodinger 
equation describing the motion of the atoms, which is proportional to the local density of 
atoms and which gives rise to the celebrated Gross-Pitaevskii equation 39|. The strength 
of this nonlinear contribution can be controlled by the scale of the confinement in the 
transverse (vertical) spatial direction. We shall make, in the following, the simplifying 
assumption that this nonlinearity is present only within the scattering region and vanishes 
within the waveguides. We furthermore assume that the waveguides are perfectly uniform, 
and that the two-dimensional scattering geometry can be described by perfect "billiard" 
potentials which combine a vanishing potential background within the waveguides and the 
scattering region with infinitely high hard walls along their boundaries. These assumptions 
considerably simplify the analytical and numerical treatment of the problem, and allow for 
the identification of well-defined asymptotic scattering states within the waveguides. Two 
such billiard configurations are shown in Fig. [1] 

The dynamics of this matter-wave scattering process is then well modeled by an inho- 
mogeneous two-dimensional Gross-Pitaevskii equation j4o| 



2m \ i 



Mr] 



+ V{r)+g^\^{r,t)\' 



^{r,t) + Sir,t) (1) 



with r = {x,y). Here m is the mass of the atoms and V{r) represents the confinement 
potential that defines the waveguides and the scattering region. The effective vector potential 
A(r) vanishes within the waveguides. Within the scattering billiard we choose it as 



A(r) = ^Be^ x (r - fq) 



^B[{x - xo)ey -{y- yQ)e^ 



(2) 



where tq = [xq, y^) represents an arbitrarily chosen reference point and e^^, e^, are the 
unit vectors in our spatial coordinate system. In the presence of an harmonic transverse 
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Figure 1: Shapes of the biUiards a and b under consideration, plotted together with the density of a stationary 
scattering state. We indicate, in addition, the widths W and W of the incident and the transmitted 
waveguides, respectively, as well as the horizontal position at which the incident and reflected parts of 
the scattering wavefunction are decomposed in transverse eigenmodes of the waveguide. The semiclassical 
average that is undertaken in order to obtain the mean retro-reflection probability involves an average of 
the retro-reflection probability within a finite window of chemical potentials /i for different incident channels 
and different locations of the circular and the lower semicircular obstacle in the case of billiards a and b, 
respectively. 



(vertical) confinement with oscillation frequency u± = uj±{r), the effective two-dimensional 
interaction strength is given by g{r) = 4:\/27ras/a±{r) with a±{r) = ^yh/[muJ±{r)] where 
denotes the s-wave scattering length of the atoms. As stated above, we assume that g is 
constant within the billiard and vanishes in the waveguides. 

The source amplitude S{r, t) describes the coherent injection of atoms from the Bose- 
Einstein condensate within the reservoir trap. Assuming that only one transverse eigenmode 
in the waveguide is populated, we may write S as 

S{r,t) = SoXi{y)S{x - Xi)exp (3) 

where Xiiv) denotes the normalized wavefunction associated with the transverse eigenmode 
with the excitation index i, characterized by the energy Ei, into which the source injects the 
atoms from the condensate (typically one would attempt to achieve coherent injection into 
the transverse ground mode, with i = 1, in an atom-laser experiment j^). xl represents an 
arbitrary longitudinal coordinate within the waveguide (which, without loss of generality, 
is assumed to be oriented along the x axis) and fi is the chemical potential with which the 
atoms are injected into the waveguide. Making the ansatz 

^(r,t) =^(r,t)exp ( -^/it ) (4) 



we obtain 



ih^^^ir,t) = (H-fi)^{r,t)+g^\^ir,t)\^^p{r,t) + SoXiiy)5ix-XL) (5) 



with the single-particle Hamiltonian 



H 



1 
2m 



'h 



V-A(r) 



+ V{r) 



(6) 



The time evolution of the scattering wavefunction can be considered to take place in 
the presence of an adiabatically slow increase of the source amplitude 5*0 from zero to a 
given maximal value. In the absence of interaction, this process would necessarily lead to a 
stationary scattering state, whose decomposition into the transverse eigenmodes within the 
waveguides allows one to determine the associated channel-resolved reflection and transmis- 
sion amplitudes. In the special case of a perfectly uniform waveguide without any scattering 
potential and in the absence of the vector potential A, this stationary state is given by 



hp\[ii) 



(7) 



where p\{^Ji) = V2m/i — Ei denotes the longitudinal component of the momentum associated 
with the transverse mode Xi- Such a stationary scattering state is, in general, not obtained 
in the presence of interaction. Indeed, a finite nonlinearity strength g may, in combination 
with a weak scattering potential, lead to a permanently time-dependent, turbulent-like flow 



pel 

across the scattering region [41|, |42|, |43|, |27|, |40| , which in dimensionally restricted waveguide 



geometries should correspond to a loss of coherence on a microscopic level of the many-body 
scattering problem 40 . 

In the following, we shall restrict ourselves to rather small nonlinearities for which we 
still obtain, in most cases, stable quasi- stationary scattering states within the billiard under 
consideration |44|]. In the subsequent two sections, we shall develop a semiclassical theory 
for the self-consistent scattering state that is obtained as a solution of Eq. (|5]). Section |3] 
focuses on contributions related to coherent backscattering, while loop corrections in next- 
to leading order in the inverse number of energetically accessible channels are taken into 
account in section HI 



3. Semiclassical theory of nonlinear coherent backscattering 

3.1. Coherent backscattering in the linear case 

The key ingredient of a semiclassical description of this nonlinear scattering process is 
the representation of the retarded quantum Green function 

G,{vy,E)^{v\{E-H, + iQ)-^\v') (8) 

in terms of all classical (single-particle) trajectories (p^, q^)(t) within the billiard, indexed 
by 7, that propagate from the initial point r' to the final point r at total energy E. Here, 
we deliberately exclude the vector potential A(r), i.e. the underlying Hamiltonian is given 
by 

Ho = ^ + V{v) (9) 
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where p = —ikV represents the quantum momentum operator. The semiclassical repre- 
sentation of the Green function can be derived from the Fourier transform of the quantum 
propagator in Feynman's path integral representation, which is evaluated in the formal limit 



h ^ using the method of stationary phase. It reads j45 



Go(r,r',E) = ^A^(r, r',E)exp 



.IT 



Here, 



S^{r,r',E) 



(10) 



:ii^ 



is the classical action integral along the trajectory 7 {T^ denotes the total propagation time 
from r' to r), represents the integer Maslov index that counts the number of conjugate 
points along the trajectory (which, in a billiard, also involves twice the number of bouncings 
at the walls, in addition to the number of conjugate points inside the billiard), and 



A(r,r',E) 



271 



detD^S^{r,r',E)\ 



\/2TTih 

is an amplitude that smoothly depends on r and r', with 

9(p',r',T) 



(12) 



detD^S^{r,r',E)\ 



det 



d{r,r',E) 



(13) 



the Jacobian of the transformation from the initial phase space variables (p', r') and the 
propagation time T to the final and initial positions (r, r') and the energy E. 

The presence of a weak effective magnetic field is now incorporated in a perturbative 
manner using the eikonal approximation. As shown in Appendix B, this yields the well- 
known modification of the Green function 



G(r, r', E) = Y, r', E) exp | [S,{r, r', E) - 0,(r, r', E)] - z^/i. 



(14) 



with (f)j{T, r', E) = — 9?^(r, r', E) — (p^{r, r', E) and 



cp^{r,r',E) = - / p.,{t) ■ A[q^{t)]dt . 



m 



(15) 
(16) 



where the integration is peformed along the unperturbed trajectory q^(t). While the latter 
(diamagnetic) contribution tp^^'^ gives only rise to a spatial modulation of the effective po- 
tential background within the billiard, the former (paramagnetic) contribution ip^ explicitly 
breaks the time-reversal symmetry of the system and plays a crucial role for the intensity of 
coherent backscattering. 
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This expression for the Green function can be directly used in order to construct the 
scattering state ip{r) that arises as a stationary solution of Eq. ([5]). We obtain 



ip{r) = So / G[r,{xL,y'),fj]Xiiy')dy' 



(17) 



where Xi represents the energetically lowest transverse eigenmode within the waveguide. 
Assuming billiard-like waveguides with a vanishing potential background and infinitely high 
hard walls along their boundaries, the nth normalized transverse eigenmode [n > 0) is given 
by 



Xn{y) 




2 . fPny 



i /A 



exp ( -Pny 



exp ( -^Pny 



ioTO<y<W (18) 



and Xn{y) = otherwise, = rnih/W is the quantized transverse momentum and W 
represents the width of the waveguide. We can therefore write 



^(r) 



Sq / Tlh 



{G [r, {xL,Pi),fj] - G [r, (xl, -pi), /i] } 



where 



G[r, {x\p'),E] 



w 



-^y^" G[r, (x',i/'),^]exp [^jp'^y' ]dy' 



(19) 



(20) 



denotes a partial Fourier transform of G[r, (x', y'), E\. 

Inserting the semiclassical expression ([Hj) for the Green function G, this partial Fourier 
transform can again be evaluated using the stationary phase approximation. The stationary 
phase condition yields (piy)j^[r, (x', y'), -E] = p^, i.e. p'y should be the y-component of the 
initial momentum of the trajectory. The integration over y' yields the prefactor ^/^jrihja 
with 

-,S^[v,[x ,y),E] = ..A 7:^1 • (21) 



a = 



dy 
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d[v, {x',y'),E] 



Combining it with the prefactor y | detD'^S^\ according to the expression f ll3l) and with the 
other prefactors that are contained within the amplitude A^, we finally obtain 



G(r, z', E) = Y, A(r, z', E) exp \ \S,{v, z', E) - 0^(r, z', E 



2' 



(22) 



with 



A(r,z',E) 
:S,(r,z',E) 



27r\/i 



9[p',(x',2/'),T] 



\/2-Kih 



det 



9[r,(x',p'),E] 



l:^5^(r,r',E)<0 
: otherwise ' 

S, {r, [x', y'^ir, z', E)] , E} + ^^(r, z', i?) 
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(23) 

(24) 
(25) 



and 4)^{r,z',E) defined according to Eq. f lB.9|) . where the initial phase-space point of the 
trajectories 7 is given by the combination z' = {x',Py) and y'^ir, z', E) denotes the resulting 
initial y coordinate. 

Channel-resolved reflection and transmission amplitudes can now be computed by pro- 
jecting ip onto the transverse eigenmodes of the waveguides. This involves again a partial 
Fourier transform of the Green function, this time in the final coordinate. In particular, the 
reflection amplitude into channel n is obtained from 



= 5, 



w 

Xn{y)'^{xL,y)dy 

-G[{Xl, -Pn), {xL,Pi),^A + G[{Xl, -Pn), (Xl, -Pi) , i4 

7ih 



'W L 



G{z„, , z+, fi) - G{z+, z+, fi) - G{z^ ,z^,fi) + G(z+, z^ , fi) 



with 



G[ix,Py),r',E] = 
d[ix,Py),z',E] ^ 



w 



y2nhJo 
1 



G[{x,y),r',E]exp ( -^Pyy ]dy , 



V2TThJo 



w 



G[{x,y),z',E]exp ( -^Pyy ]dy 



where we define 



z±^ 



{xl, ±Pn) for incoming trajectories (with p'^ > 0) 
{xl, TPn) for outgoing trajectories (with Px < 0) 



Similarly as for G, the semiclassical evaluation of this Fourier transform using Eq. 



G{z, z', E) = J2 z', E) exp \ \s^{z, z', E) - 0^(z, z', E) 



IT- 



(26) 



(27) 

(28) 
(29) 

(30) 
yields 

(31) 



with 



A^{z,z',E) 
5'-y(z, z', E) 



2m 



dM,p'y),{x',y'),T] 



det 



d[{x,Py),{x',p'y),E] 



\/2TTih 

S-f {[x,y^{z,z',E)],z',E} -pyy^{z,z,E) 



l:^S,ir,z',E)<0 
: otherwise 



(32) 

(33) 
(34) 



and (j)^{z,z',E) defined according to Eq. fIB.Qp . where the final phase-space point of the 
trajectories 7 is given by the combination z = {x,py) (and y^ is the final y coordinate). 
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From Eq. (13T|) it becomes obvious that subtle interferences between different classical 
trajectories may give rise to channel-resolved reflection and transmission probabilities that 
strongly fluctuate under variation of the incident chemical potential fi = E. Those fluc- 
tuations generally cancel, however, when performing an average within a finite window of 
chemical potentials. Specifically, the calculation of P involves sums over pairs of tra- 
jectories 7 and 7', whose contributions contain phase factors that depend on the difference 

— Sy of the associated action integrals. These differences strongly vary with the chemical 
potential fi unless the two trajectories 7 and 7' are somehow correlated. 

An obvious correlation arises if the two trajectories happen to be identical, in which case 
the phase factor is unity. In the framework of the diagonal approximation, we only take into 
account this specific case, i.e., we approximate the double sum J^yy by a single sum 
where 7' is taken to be identical to 7. The energy average (|'?/'nP) of \4'n\'^ is then given by 



Sn 



W 



2 r 



G'(z+,z+,/i) 



+ 



+ 



(35) 



(36) 



G(z,z',E) 



V ( A,{z,z\E) 



with ^ 

As shown in Append ix C[ this sum is evaluated using the generalized Hannay-Ozorio 
de Almeida sum rule [48, 49]. Defining by td the "dwell time" of the system, i.e. the mean 
evolution time that a classical trajectory spends within the billiard before escaping to one 
of the waveguides, and introducing the "Heisenberg time" as = mQ/h where Q denotes 
the area of the billiard, we obtain [see Eq. flC.lSp ] 



G(z,z',E) 



mW 



1 



2nh^J th ^2mE - p'^ ^2mE - p'^ ' 
Inserting this expression into Eq. fl36|) and defining 

pI{E) = ^2mE - pI = ^/2mE - {mrh/Wy 



(37) 



(3^ 



as the longitudinal component of the momentum that is associated with the transverse mode 
Xn finally yields 



(l^nl')d 



mSo 



rnpMplif^) 



(39) 



This expression can be used in order to determine the steady current jn of atoms that are 
reflected into channel n, according to 



Jn 



m 
10 



(40) 



Dividing it by the incident current which is derived from Eq. as 



,_ m\So\^ 

we obtain the reflection probabihty into channel n as 

r-ai = jn/f = Td/th • 



(41) 



(42) 



The same reasoning can be applied to the outgoing waveguide on the other, transmitted 
side of the billiard. Again we obtain t„j = td/th as the probability for transmission into 
the transverse channel n of the outgoing waveguide, even if its width W is different from 
the width W of the incoming guide. The total reflection and transmission probabilities R 
and T are then simply related to the numbers of open channels Nc and Nc in the incoming 
and outgoing waveguide according to R = Nc^dI^u and T = Nc'Td/'Thi where we evaluate 
Nc = 2W/XdB and Nc = 2iy/AdB in the semiclassical limit, with AdB = 27ih/ \/2'mfj, the de 



Broglie wavelength of the atoms. We can furthermore use the general expression [29|, |30 



{W + W)v 



(43) 



for the mean survival time of a classical particle propagating with velocity w in a chaotic 
billiard with area Q that contains two openings of width W and W, which yields 



To AdB/2 



th W + W Nc + Nc 



(44) 



We then arrive at the intuitive results R = W/{W + W) and T = W/{W + W), i.e. the 
total reflection and transmission probabilities are simply given by the relative widths of the 
corresponding waveguides. 

The diagonal approximation therefore yields predictions for reflection and transmission 
that are expected for incoherent, classical particles in a chaotic cavity. It represents in 
leading order in the inverse total channel number {Nc + Nc)~^ the contributions for all 
channels on the transmitted side, and for all reflected channels except for the channel n = i 
in which the matter-wave beam is injected into the billiard. In this incident channel, there 
is another, equally important possibility to pair the trajectories 7 and 7' in the double 
sums that are involved in the calculation of {ipil"^: 7' can be chosen to be the time-reversed 
counterpart of 7, the existence of which is guaranteed by the time-reversal symmetry of Hq. 

Consequently, Eq. (^5^ has to be corrected for the special case n = i according to 



(45) 



where the "crossed" or "Cooperon"-type contribution 



2 r 



+ 



^(z. ,z. ,/i) 



+ ( G (z+, z. , /i) G(z. , z+, /i) ) +{G (z,^ , z+, /i) G(z+, z- , fx) 
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(46) 



contains all those combinations of trajectories for which 7' is the time-reversed counterpart of 

7. Obviously, the action integrals and Maslov indices JI^ are identical for the trajectories 
7 and their time-reversed counterparts. This is not the case, however, for the modification 

(p^ of the action integral that is induced by the vector potential, whose paramagnetic part 
Ip^ [Eq. f lT^ ] changes sign when integrating along the trajectory 7 in the opposite direction. 
We therefore obtain 



g\z',z,E)G{z,z',E)) =J2 



A^{z, z', E) 



exp 



(/?^(z,z',E) 



(47) 



To provide some physical insight into the role of this additional phase factor, we use the 
representation ([2j) of the vector potential within the billiard. Using p-Y(t) = mq^(t) along 
trajectories 7 generated by iJo, the paramagnetic contribution to the effective action integral 
reads then 



(z,z',i?) 



B 



[q^(t) - ro] X q^(t) dt 



(4J 



2 .0 

where ro is an arbitrarily chosen reference point. Within the billiard, the trajectories 

(p^, q^)(t) can be decomposed into segments of straight lines that connect subsequent re- 

(1) AN-i) 

.7 5 • • • 5 47 

r, the initial and final points of the trajectory, we rewrite 



flection points at the billiard boundary. Denoting those reflection points by q' 
and defining q^'^^ = r' and q^^^ ~ ^ ^^^^f^oi ^-t^A r^^^^^fc fi^^ fvoi^^f^ 



Eq. 



where 



as 



N 



ip^{z,z',E) = Bj2(^j=BA 



1 

- 

2 



0-1) 

7 



ro) X 



(49) 



(50) 



is the directed area of the triangle spanned by the reflection points q^ and q^ '' as well as by 
the reference point rg. Quite obviously, ip^ is independent of the particular choice of tq, or of 
any other gauge transformation A 1— )■ A + Vx that vanishes within the waveguide, provided 
the initial and final points r' and r of the trajectory 7 are identical or, less restrictively, lie 
both within the same, incident waveguide where the vector potential vanishes (in which case 
a straight-line integration J A ■ dq from r' to r would formally close the trajectory without 
adding any further contribution to ^^). 

The central limit theorem is now applied in order to obtain the probability distribution 
P{T-y,A) for accumulating the area A after the propagation time |29|, |30|, |46|, |47|. We 
have 

^(^7, ^) = . i,. ^ exp ) (51) 



where Vl is the area of the billiard, v is the velocity of the particle, and r/ is a dimensionless 
scaling parameter that characterizes the geometry of the system and that can be numerically 



computed from the classical dynamics within the billiard as described in Appendix D This 
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distribution is now used to obtain an average value of the magnetic phase factor according 
to 

2?. 

' ' ~ ' ' ' (52) 



exp 



J P{T^, A)exp (^BAj dA = exp (^-^ 



with 



(53) 



the characteristic time scale for magnetic dephasing. 

With this information, we can now follow the derivation of the Hannay-Ozorio de Almeida 



sum rule, as explicated in Appendix C in order to evaluate the expression P7|) . with the 
only complication that each contribution in the sum over trajectories needs to be weighted 
by the "dephasing" factor exp(— T^/r^). This yields 



G (z',z,E)G'(z,z',E) 



Hence, we obtain 



+ 



/2mE - pI ^2mE - p'^ 



mSo 


2 






hp\{fi) 









in very close analogy with Eq. (139|) . which altogether yields 

uiSq 



hp\{fi) 



1 



1 + td/tbJ Th 



(54) 



(55) 



(56) 



This gives rise to an enhanced probability for retro-refiection into the incident channel n = i, 
namely 

1 



1 + 



1 



with 



1 + td/tb 



1 



1 



(57) 



(5^ 



^/ryvTo^W^ 

as compared to reflection into different channels described by Eq. fH21) . which is the char- 
acteristic signature of coherent backscattering. Note that, due to conservation of the total 
flux, increased retro-refiection for n = i implies decreased reflection or transmission into 
other channels n ^ i. This will be subject of Section H] below. 

The above prediction (1571) is expected to be valid for chaotic cavities in the semiclassical 
limit of small h (i.e. of a small de Broglie wavelength as compared to the size of the scattering 
region) and in the limit of small widths of the leads. Leads of finite widths, as the ones that 
are considered in the scattering geometries shown in Fig. [1], will give rise to non-universal 
corrections to Eq. fl57|) that are related to short reflected or transmitted paths. In particular, 
the presence of self-retracing trajectories, which are identical to their time-reversed coun- 
terparts, affects the probability for retro-reflection due to coherent backscattering, as those 
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trajectories are evidently doubly counted in the addition of ladder and crossed contributions. 
Hence, the enhancement of this retro- reflection probability with respect to the incoherent 
ladder background (H2l) will, in practice, be reduced as compared to Eq. (157|) . due to the 
presence of short and therefore semiclassically relevant self-retracing trajectories. 

3.2. Diagrammatic representation of nonlinear scattering states 

We now consider the presence of a weak interaction strength g > Oin the Gross-Pitaevskii 
equation ([5]). As a consequence, the scattering process becomes nonlinear and the final 
(stationary or time- dependent) scattering state may depend on the "history" of the process, 
i.e. on the initial matter-wave population within the scattering region as well as on the 
specific ramping process of the source amplitude. We shall assume that the scattering 
region is initially empty (i.e., V'(r, i) = for t — )• — oo) and that the source amplitude 5*0 is 
adiabatically ramped from zero to a given maximal value 5*0, on a time scale that is much 
larger than any other relevant time scale of the scattering system. This adiabatic ramping 
is formally expressed as S'o(t) = Sofit/tji) where /(r) is a real dimensionless function that 
monotonically increases from (for t — i- — oo) to 1 (for t — )■ oo) and t/j — )• oo is a very large 
ramping time scale. Redefining ^/'(r, t) = f{t/tR)tlj{r,t) and neglecting terms of the order of 
l/tfi, we obtain from Eq. ([5]) 

ih^^H^, t) = {H- /x)V^(r, t) + Sir, t) (59) 

as effective Gross-Pitaevskii equation for with 

^(r, t) = SoX.iy)5ix - xl) + t)\^{v, t) (60) 

and g{t) = P{t/tR)g- For weak enough nonlinearities g and long enough ramping time 
scales tR, Eq. f l59|) can be considered as describing an effectively linear scattering problem 
the source term of which is gradually adapted according to Eq. f l60|) . We can therefore 
express the time-dependent scattering wavefunction as 

^(r,t) = J d^r'G{r,r',fi)S{r',t) (61) 

where G = {fi—H+iO)~^ is the Green function of the linear scattering problem [see Eq. f|T^ ]. 
In the limit of long evolution times t — oo, we thereby obtain 

^{r) = SoJ G[r,{xL,y'),fiMy'W + J rfVG(r, r', /i)^?|^|^(r')PV^(r') (62) 

as self-consistent equation for the scattering wavefunction, which generalizes the expression 
ffT7|) obtained for the linear case. 

In rather close analogy with the numerical procedure that is employed for computing a 
stationary scattering state, we can construct a self-consistent solution of Eq. (162|) by starting 
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with the expression (|T7|) for the hnear case and by iteratively inserting the subsequent 
expressions obtained for ip{r) on the right-hand side of Eq. (I62|) . This naturally leads to a 
power series in the nonlinearity, 



(63) 



n=l 



where tp^^^ (r) represents the solution of Eq. f|T7|) . i.e. the scattering state of the noninteracting 
system. 

It is instructive to evaluate the semiclassical representation of the first-order correction 
to the linear scattering wavefunction ip^'^\ given by 

5^(1) (r) = — [ dVG(r,r',/i)|^(°)(r')|'^(°nr') ■ (64) 
2m J 

Using the expression (|T9|) for the scattering state of the noninteracting system, we obtain 



dip 



h"^ So h^h 



E 



Z/1Z/2Z/3 



X / c^VG(r,r^/i)G(r',zr^/x)G'(r',zr^/x)G(r',zr^/x) 



(65) 



with zf^ = as defined in Eq. (!30l) . Inserting the semiclassical expansion for the Green 
function, given by Eqs. (HM and yields 



fi? Sq /vr/i 

2^TV ^ 



Z/1Z/2Z/3 



j d'r'Y, E Ao(r,r',/i)A,(r',zr,/i)A;(r',z-,/x)/l,3(r',zr,/x) 



70 71 1 72, 73 



X exp |- [5^0 (r, r', /i) + ^^^(r', z^^^, /i) - ^^^(r', zf , /x) + ^^3(r', zf , /i 
X exp (-^ [0,„(r,r',/x) + 0^^(r',zr,/x) - 0^^(r',zr,/i) + 0,3(r',zf ,/x)] 
X exp 



tn / \ 

-y (/^70+/^7l-/^72+/^73) 



(66) 



where the indices 70 and 7^ (£ = 1, 2, 3) represent trajectories that connect r' and r as well 
as zj'* and r', respectively. 

Neglecting, as done in Section |3l the modification of the trajectories 70 and 71/2/3 due to 
the presence of the weak magnetic field, a stationary-phase evaluation of the spatial integral 
in Eq. flM]) yields the condition 



P70 ^'^ + P72 (j"'' ' = P71 (^'^ > /^) + P73 (^'' ' /^) • 



(67) 
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Noting that all involved momenta are evaluated at the same spatial point r', this condition 
is satisfied if and only if 



or 



or 



P7o(r,r',/i) = p^^(r',zf and p\^{r' , z"^'' , fi) = pi^,^{r' , z^^'', fi) (68) 
P70 (r , r , ^) = (r', z'^' , fi) and p^^ (r', , /i) = p^ ^ (r, z^^ , fi) (69) 



p;(r,r',/i) = -p;(r',zr,/i) and p;^(r', zf, ^) = -p;(r', zf , ^) (70) 

holds true. The cases (1^ and (jSS]) are essentially equivalent and imply , in case (^E\i [or in 
case fl69|l ]. that the trajectories 72 and 73 (or 72 and 71) are identical and that 70 represents 
the direct continuation of the trajectory 71 (or 73) from r' to r. This latter condition 
determines the stationary points of r', which have to lie along the trajectories from z^^ (or 
zf ) to r. 

Case fl7UI) is more involved. It implies, on the one hand, that the time-reversed counter- 
part of trajectory 73 represent the direct continuation of trajectory 71 (using the fact that 
the scattering system under consideration is, in the absence of the magnetic field, invariant 
with respect to time reversal), which determines the stationary points of r' along refiected 
trajectories from z^^ to z^^. On the other hand, 70 represents a part of the time-reversed 
counterpart of trajectory 72, which necessarily implies that the point of observation r has 
to lie along 72. This latter condition generally represents an additional restriction of the set 
of stationary points in Eq. fl66|) (namely that r' lie on the continuation of a trajectory from 
zj'^ to r), which substantially reduces the weight of contributions resulting from case flTOjl 
as compared to those emanating from cases fl^Sl) and fIBIJl) . An exception to this rule arises 
if the point of observation r is identical with or lies rather close to z^^, in which case all 
contributions resulting from Eqs. flBSjl - flTO]) are of comparable order. 

In full generality, we can express the first-order correction to the linear scattering wave- 
function in the semiclassical regime as 

(r) = 2(5V^f ^ (r) + 6^^^^ (r) (71) 

where 6ilj[^\r) and 6ilji^\r) contain the contributions that respectively emanate from the 
cases (1681) . fl69p as well as from the case fl70l) . Considering an observation point r that lies 
deep inside the billiard, we neglect ^/'c^^(r) for the moment. The expression for 6tlj^^\r) can 
be cast in a form that is, apart from a source- dependent prefactor, exactly equivalent to 
the first-order term in the Born series of a perturbed Green function, where the effective 
perturbation Hamiltonian SH corresponds here to the density |'?/'''°^(r)|^ of the noninteracting 
scattering wavefunction as evaluated by the diagonal approximation, i.e. to 



J^|A^(r,z+,^)|V^|A^(r,z, ,/i)'^ 



(72) 
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In close analogy with the first-order modification (]B.6|) of the semiclassical Green function 
in the presence of a weak perturbation, we then obtain 




xA^{r, zj', /i) exp <j ^ ^^(r, zj', /x) - 0^(r, z^, /x)] - 



(73) 



Sip^l^ and (^-i/^c^^ shall, in the following, be termed "ladder" and "crossed" contributions, 
respectively. 

To illustrate this point, it is useful to introduce a diagrammatic representation for this 
nonlinear scattering problem. Following Ref. 3J|, we represent by ^ and ^ 



the Green function G(r, r', /x) and its complex conjugate G*{r, r', /x), respectively. The (four- 
legged) vertex H represents a scattering event of ip at its own density modulations, described 
by the second term of the right-hand side of Eq. fl62l) . and I denotes the corresponding vertex 
for i/j*, appearing in the complex conjugate counterpart of Eq. fl^ . The source is depicted 

by the vertical bar I, i.e. I ► represents the scattering wavefunction of the noninteracting 

system, given by the convolution of the Green function with the source. We can then express 
Eq. (162|) and its complex conjugate as 



(74) 




(75) 



where ^"^^^ and ^"^^^ respectively represent the self-consistent stationary scattering 
wavefunction ip{r) of the nonlinear system and its complex conjugate '?/'*(r). Going up to the 
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second order in the power-series expansion (1631) . we obtain the diagrammatic representation 




+ 0{g') . (76) 



The semiclassical evaluation of the first-order term according to Eqs. fl7T]) and 
neglecting the contribution of 5il)^\ can be expressed as 

^ ^ 2 I— ►H—^ (77) 

I It 

in diagrammatic terms. In close analogy with the corresponding ladder diagrams in disor- 



dered systems 32|, |33|, |3J| , the parallel arrows ^ symbolize the semiclassical evaluation 
of G*G in the diagonal approximation, with G and G* following the same trajectories that 
connect a given initial with a given final point. The diagram ' ►H ►, on the other hand, 
indicates that the nonlinearity event takes place along a continuous trajectory that connects 
the source with a given final point at the end of the arrow. As already discussed above, the 
factor 2 in Eqs. f lTT]) and fl77|) originates from the two equivalent conditions f pS]) and flB^ . 
In other words, the red arrow on the left-hand side of Eq. (1771) can be paired with either one 
of the two incoming black arrows. 

3.3. Ladder contributions 

It is suggestive to pursue the analogy with the Born series of a linear Green function 
and to introduce a modified Green function G^ (the i stands for "ladder contributions"), 
symbolized by in which the contribution of the density-induced perturbation is 

summed up to all orders in the nonlinearity g. The Dyson equation that this Green function 
satisfies is represented as 



>+2 ► +4 ^ +... 

tt t It 

> + 2 ^ . (7J 



ir 
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Applying the stationary phase approximation, the explicit expression for this modified Green 
function reads, in analogy with Eq. f lB.8|) . 



Gi{r,r',fi) = ^ A^(r,r',/i)exp |^ [^^(r, r', /i) - 0^(r,r',/i) - x-y{r,r' , fi)] -i^Hj^ 



(79) 



with x^(r, r',/i) = 2g{h^/2m) J^"' \ilj^^'^[q^^{t)]\'^dt. On this level, the nonlinearity therefore 
induces an effective modification of the action integral along the trajectory 7, in close analogy 
with the change in action for the dynamics in the presence of a weak static disorder potential 



50|. This modification, however, does not at all affect the calculation of mean densities 
within the billiard using the diagonal approximation: evaluating the wavefunction ^p{r) 
according to Eq. ( [T9|) with G being replaced by Ge, we would essentially obtain |'?/'(r)|^ = 
iV'^^Hr)!^, the latter being given by Eq. (172|) where the phases x-y appearing in Eq. (1791) drop 
out. 

The same reasoning applies if we replace ^^^^ by ip in the definition of the nonlinearity- 
induced modification of the effective action associated with the trajectory 7, i.e., we (re-)de- 
fine ^ ^ 

X^(r,r',/i) = (7- / ' mq-yimldt (80) 
Jo 

and use this expression in the definition of Ge according to Eq. (179|) . This amounts to 
replacing the diagrammatic representation (ITHj) by 



>+2 



tt 



which, when being expanded in powers of g and evaluated using the stationary phase ap- 
proximation, involves all possible ladder- type (parallel) pairings of G and G*, i.e.. 



► = ^+2 ^ +4 

It tt tt 

+4 ^ +4 ^ +0{g') (82) 

up to second order in g. The mean density within the billiard as evaluated using the diagonal 
approximation is then given by 



|2 _ IC |2^^ 

Id - l-^ol 



^|A^(r,z+,/i)f + ^|A^(r,z. ,/i)| 
. 7 7 
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(83) 



as in the case of the hnear scattering problem [see Eq. (172|) ] 51 



It is worthwhile to calculate the energy average of the density within the billiard using 
the Hannay-Ozorio de Almeida sum rule 



49|. As shown in Appendix C we have [see 



Eq. fini3|) 



E 



Ay{r,z',iJ,)\ 



rn?W td 



This eventually yields 



(I^WI'>d 



mSn 



h 



■r 



i5) 



when being expressed in terms of the incident current = m\ S ol"^ /[h'^p\{fi)]. The mean 
density is therefore obtained from an equidistribution of the population in the case of a 
stationary flow, which is given by the ratio of the feeding rate j' and the decay rate r^^. 

3.4- Crossed contributions 

As seen above, the nonlinear ladder contributions vanish on average. However, we have 
so far neglected the influence of terms arising from the association of trajectories according 
to the remaining (and less intuitive) case f lTU]) . As was argued above, the contributions 
of such terms to the local density is generally suppressed with respect to the ladder-type 
contributions arising from the cases fl^ and (En]), due to the fact that case ( 170]) requires 
not only the time-reversed counterpart of trajectory 73 to represent the direct continuation 
of trajectory 71, but also that the point of observation r lie on the trajectory 72 connecting 
the source with the interaction point r'. In the case of retro-refiection into the incident 
channel, however, where r lies directly at the location of the source, this latter condition is 
satisfied by default, and we should therefore expect a finite contribution from this "crossed" 
association of trajectories to the probability of coherent backscattering. 

It is instructive to first compute the influence of such crossed terms in linear order in 
the nonlinearity. We evaluate for this purpose the remaining term 6tlji^\r) in Eq. fl7T]) that 
is associated with the case f lTU]) . The requirement that the time-reversed counterpart of 73 
represent the direct continuation of 71 allows one to apply the stationary phase approxima- 
tion in order to evaluate the spatial integral in Eq. ( 166]) . In close analogy with Eq. ( 173]) . we 
then obtain a single sum over all trajectories 7 that connect the initial phase-space point 
z^^ with the final point zj'^ {1^1,1^3 = ±1) both being associated with the incident channel 
Xi{y)- An important extension as compared to the structure of Eq. (I73p is provided by the 
paramagnetic contribution (ITB]) to the effective action integral, which changes its sign under 
the time-reversal of the trajectory 73. 

Calculating the overlap of 6ilji^\r) with the incident channel, we obtain the associated 
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first-order correction to the backscattering amplitude as 



W 



Sr 



Xi{y)Hc\xL,y)dy 



(86) 



/i) 



X exp 



^[5,(zr^z^,/.)-0^(zr^zr^/i) 



TT- 



where we define 



2^^^ / C(°)[q,(t)]exp^--^,[zr,q,(t),/i]|.rft 



2z. 



^7) 



72 



2^- r + ^' 



72 



. _ . 



as "crossed density" within the bilhard. The latter quantity can be interpreted as the 
semiclassical evaluation of 



l^ol 



,TTh 

w 



G*{r, z+ /i)G(z+, r, /i) + G*(r, z- , /i)G(z. , r, fi) 



(89) 



within the diagonal approximation. In contrast to the actual density within the billiard, 
given in leading order by the expression ( 183|) , is, in general, not invariant under gauge 

transformations A 1— )• A + Vx of the effective vector potential A(r), due to the presence of 
the phase factors containing the paramagnetic contribution to the effective action integral. 
The combination of those phase factors with the corresponding one arising in Eq. fl87p . 
however, gives rise to an overall expression that is invariant under gauge transformations. 

To verify this, we introduce for each point r within the billiard a straight-line trajectory, 
denoted by the index u, that connects this point to a fixed reference point within the 
incident lead, given, e.g., by = {xl, W/2). This straight-line trajectory can be defined as 



r+—{rL 



(90) 
(91) 



with T^^ = m|r2, — r\j ^2m[i |5j]. We now define for each "incident" trajectory 7, i.e. which 
connects a phase-space point z within the incident lead to a spatial point r within the 
billiard, its "completion" as 7 = a; o 7. In physical terms, 7 traces the motion of a particle 
that follows 7 and is then scattered back to the incident lead due to the presence of a local 
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perturbation within the bilhard (a point scatterer) at position r. We obviously have the 
relation 

^^(rL, r, z, /i) = ifu^iTL, r, /i) + ^^(r, z, /i) (92) 

for the paramagnetic action integral along the trajectory 7. As integrations of p(t) ■ A[q(t)] 
along paths that are entirely contained within the incident lead obviously vanish due to 
the local absence of the vector potential, we can state that Tp;y{rL, r, z, fi) is invariant under 
gauge transformations. Analogously, a trajectory 7' that leads from a spatial point r within 
the billiard to a phase-space point z within the incident lead is "completed" as 7' = 7' o a; 
with the associated paramagnetic action integral 



(93) 

The crossed density fl88|) can therefore be re-expressed in terms of such completed tra- 
jectories 72 through 



c^^^ (r) exp 



(94) 



where its gauge-invariant part is introduced as 



,z+,/i) exp 



72 



-- ^<^^2(rL,r,z+, 



/i) 



rL,r,z. ,/i)| exp 



72 



. _ ; 

--^V^^2(rL,r,z. ,/i) 



(95) 



As in the case of "ordinary" backscattering trajectories, the energy average of the param- 
agnetic phase factor of 7 yields, in analogy with Eq. fl5^ . 



exp 



2z_ ■ 
-- ^V^^(rL,r,z,/i) 



exp I 



(96) 



with tb the characteristic time scale associated with the magnetic field, defined by Eq. (153|) . 
We neglect in this expression the contribution of to the total propagation time of 7 (which 
is, in fact, canceled in the nonlinear diagrams contributing to the backscattering probability 
to be discussed below, as the latter involve, by construction, flux integrals along closed 
paths) and assume T^^ ~ Ty. In perfect analogy with the derivation of the energy- averaged 
density within the billiard, we then obtain [see Eqs. fIC.Sp and flC.lSp 



(cW(r)> 



h 



m 



(97) 



for the energy average of the gauge-invariant part of the crossed density. 

This expression can be used in order to evaluate the first-order correction to the crossed 
contribution of the nonlinear backscattering probability according to 
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(98) 



where P)c*^^ represents the hnear crossed contribution as defined in Eq. (H6l) . Within the 
diagonal approximation, we obtain 



h2m 



2i 



A(zr,<s/i) 



rft(c(°)[q,(t)]> 



X (exp <j -— ¥?^[z,^%q^(t),rL,/i] }> + exp <j — y^^^fr^, q^(t), z^S /i] \ )(99) 



2i 



where, in the second fine of Eq. fl99|l . we account for the fact that each trajectory 7 can 
be paired with itself as well as with its time-reversed counterpart, the latter giving rise to 
a different paramagnetic phase factor. For both possibilities of the pairing, the remaining 
piece of the trajectory 7, respectively connecting q^(t) with z^^ as well as z^^ with q^(t), can 
be "completed" by combining it with the straight-line trajectory uj from q'y(t) to that is 
introduced through the factorization (l94l) . 

We can now perform the energy average of the paramagnetic phase factors according to 
Eq. (l96i) . taking into account that the effective propagation times of the pieces of trajectories 
under consideration equal T^ — t as well as t, respectively, for the two phase factors appearing 
in the second line of Eq. fl99p (the additional contribution of the straight-line trajectory uo to 
the total propagation time is neglected). For both phase factors, this gives rise to integrals 
that are straightforwardly evaluated as 









--) 







dt 



tb 



exp 



^b) 



We therefore obtain 



h 2m 



Tih 



) E 2-«E 



which, after applying the Hannay-Ozorio de Almeida sum rule [see Eq. 
as 



mSn 



-2 ^1 



m 



(100) 



^7 

exp I 



(101) 

Lisp ], is evaluated 



(102) 



using the expression f l97|) for the average of the crossed density (^c^^-*). As this expression is 
purely imaginary, the modification of the backscattering probability due to the presence of 
the nonlinearity vanishes in first order in ^f, as seen from Eq. fl^S]) . 

Going beyond the first order in ^f, we can express the full nonlinear coherent backscat- 
tering probability, as evaluated using the semiclassical stationary phase approximation, in 
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diagrammatic terms according to 



1 



I 1 



+rtHtt + 2 ^Ji^ + 2 M=^M^ + 0{9') . (103) 



I 



Here, ^ represents, according to Eq. (I5T]) . the modified Green function due to 

tlie inclusion of ladder contributions. All types of ladder diagrams that were discussed in 
the previous subsection 13.31 are therefore implicitly included in this representation. As in 
Eq. f l82|) . the prefactors 2 symbolize the fact that two different possibilities of pairings have 
to be counted for certain diagrams. 

In analogy with the derivation undertaken in Ref. |34y], this series of diagrams can be 
exactly summed yielding 



r 



+ 



+ 



(104) 



where we define the nonlinear crossed density 
self-consistent manner through 



and its complex conjugate 



m a 



It 



+ 2^ 



I 
I 



(105) 



(106) 
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This nonlinear crossed density can be expressed through a transport equation of the form 



C^^)(r) + 2g—\S, 



2m 



,Trh 
W 



1/1,I/2=±1 



Z/iI/2 



X J rfVc(^)(r')G,(r',r,/.)G,(r',zr,/x)G;(r,zr,/x) (107) 



which involves the modified Green function fj79|) that takes into account the average shift 
of the effective potential within the billiard due to the presence of the nonlinearity. Being 
invariant under time-reversal, the nonlinearity-induced contribution x-y fISU]) to the action 
integral does not play any role for the determination of the nonlinear crossed density. Indeed, 
applying the stationary phase and diagonal approximations in Eq. fll07p . we obtain 

C7(.)(r) = C(°)(r) + ^|5or^5^p,(r,zr,/.)f 



mW 



u=±l 



2i 



C(^)[q^(t)] exp \ --^f-ziT^, <i-/{t),fA \ dt 



:i08) 



which does not involve any reference to x-y- 

Quite obviously, the nonlinear crossed density C^^\y) is not invariant under gauge trans- 
formations of the effective vector potential. In perfect analogy with C^^\v)^ however, we can 
describe the explicit gauge dependence of C^^\y) in terms of a phase factor that contains 
the paramagnetic contribution of a straight-line trajectory u from r to the reference point 
vl within the incident lead [52|]. In analogy with Eq. fl94|) . we therefore propose 

C(9)(r) = c(^)(r)exp 
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as definition for the gauge-invariant part c^^\r) of the nonlinear crossed density, which in 
turn satisfies the gauge-invariant transport equation 
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We can now compute the energy average of c'^^\v) by assuming that it is, as the one for 
c^°-'(r) [see Eq. flW|) ]. independent of the position r within the billiard, which is to be verified 
a postenon. Using Eqs. (ESD, dHZD, f lMJ]) . (EH), and f lU:i3|l . we obtain 
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which is straightforwardly solved as 
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\ / l + ,^r^A(c(o)) 
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We are now in a position to evaluate the full nonlinear coherent backscattering prob- 
ability according to the diagrammatic representation (11041) . Denoting the linear crossed 
contribution to the backscattering probability by 
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we have, generalization of Eq. (|5 
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where we introduce 
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as contributions that result from the nonlinear diagrams in Eq. (I104p . Again, stationary 
phase and diagonal approximations are employed in order to evaluate these contributions, 
and we also use Eqs. f ll09p and f lll2p in order to express the nonlinear crossed density C-^-* (r). 
This yields for the energy average 
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(the real part of which is nonzero due to the fact that (c^^^) is complex) and 
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where we evaluate 
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Altogether, we then obtain 
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which together with Eq. ([27]) yields 
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where is the incident current. The probability for retro-reflection into the incident channel 
77, = i is then obtained as 
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This is the main result of Section [31 It essentially states that the presence of the nonlinearity 
g constitutes another dephasing mechanism in addition to the magnetic field. 

3.5. Alternative approach in terms of nonlinearity blocks 

Inspired from Refs. 32], [33], [3J], we outline, in this subsection, an alternative approach to 
determine the nonlinearity-induced modifications to the retro-reflection probability, which 
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Figure 2: Nonlinear diagram corresponding to the calculation of {ip* {Sipl'^ )) . The left-hand side shows the 
diagram together with the durations that represent the integration variables appearing in Eq. (I126p . The 
right-hand side illustrates the variable transformation leading to Eq. (|127p . In effect, the diagram can be 
cut into individual "links" the contributions of which can be determined by separate integrals. 



will become useful in the subsequent section on loop contributions. We reconsider for that 
purpose the calculation of {'^PH^ipf^)) on the basis of Eq. fl99|) . which was done using the 



expressions fl52l) and fl96l) for the average magnetic phase factors and the sum rules flC.13|l 
and flC.15|) . If we deliberately keep the occurring integrations over trajectory durations as 
they appear in the sum rules [see Eq. flC.Sp ], we obtain as an intermediate result 
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which transforms into 
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after applying the variable transformation tl = t-y, tl = Tj — t^, t'f' = T^^ that is moti- 



vated from Ref . j53| . Figure |2] illustrates this variable transformation in the corresponding 
nonlinear diagram. In effect, the transformation cuts the diagram into separate pieces of 
trajectories which we shall, as done in Ref. 53|], refer to as "links" and which are connected 
with each other at the "nonlinearity block" H. Each link gives rise to a separate integration 
yielding either td for ladder- type links t or td/{1 + td/tb) for crossed- type links -i ^ . 
Applying this rule to the diagram under consideration, we obtain again the expression fll02p 
for (^*(5^f^)). 

This reasoning can be generalized to more complicated diagrams involving more than 
one nonlinearity block. Under consideration of the sum rules (IC.lSp - (IC.lSp and of the 
combinatorial prefactors 2 arising in front of each nonlinearity block [see Eq. [H2], we can 
state the following rules: 



(1) each source I contributes a factor h/p\{jj); 
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(c) 

Figure 3: Relevant diagrams for the calculation of the nonlinear contribution to the backscattering proba- 
bility. Diagrams (a), (b), and (c) respectively correspond to the second, the third, and the fourth term on 
the right-hand side of Eq. (|104|) as well as to the second, the third, and the fourth line of Eq. (I128P . 

(2) each arrow emanating from a source I ► contributes a factor 5*0 (and each conjugate 

arrow I ► a factor Sq); 

(3) each trajectory, scattering from lead to lead or ending at a nonlinearity event within the 
bilhard, contributes a factor m'^ /{K^th); 

(4) each nonlinearity event H in the scattering wavefunction contributes a factor —igh/m 
(and each nonlinearity event H in the conjugate wavefunction contributes a factor 
igh/m); 

(5) each ladder-type link ^ contributes a factor r^; 

(6) each crossed-type link < ^ contributes a factor {I/td + 1/tb)~^- 

Using these rules, we can re-calculate the crossed contribution c[f to the retro-reflection 
intensity. In contrast to Section 13. 4[ we do not explicitly need to introduce the nonlinear 
crossed density C'-^^(r) as done in Eq. fll07p . Instead, we directly sum over all possible com- 
binations of crossed diagrams as they are depicted in Fig. |3l Together with the contribution 
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Figure 4: Basic building blocks with which all possible diagrams can be constructed. As shown in the left 
and central columns, blocks with ladder-type input can be incorporated in both the ip and the tp* branch for 
the same orientation and pairing of the trajectories outside the block region (depicted by the dashed square). 
This is not the case for blocks with crossed-type input shown in the right column, due to the mismatch of 
trajectories entering and leaving the block. Consequently, the contributions resulting from those blocks do 
not cancel each other, in contrast to the blocks with ladder-type input, but give rise to a finite modification 
of the reflection and transmission probabilities in the presence of the nonlinearity. 



of the diagram h ►!, this yields 
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which is perfectly equivalent to Eq. fll24p . 

The approach on the basis of nonlinearity blocks also provides an alternative under- 
standing why blocks with ladder-type input [i.e. where a ladder pairing is employed along 
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Figure 5: Example for a nonlinear diagram that is not accounted for in the present diagrammatic theory. 
This diagram requires the presence of two nonlinearity blocks at the same spatial location and is therefore 
not expected to be of relevance for weak nonlinearity strengths. 

the trajectory that ends at the nonhnearity event, corresponding to the cases (1^ and fl^ 
in Section I3.2[ and displayed in the left and central columns of Fig. H], do not affect the 
mean values of densities and currents of the propagating condensate. We remark for this 
purpose that the individual factor provided by each nonlinearity block is purely imaginary 
(as stated above by rule 4), with a negative imaginary part for blocks H that are incorpo- 
rated within ip and with a positive imaginary part for blocks H within ip*. Two diagrams 
that are almost identical except for the incorporation of one single nonlinearity block, which 
is placed within ip in one of the diagrams and within ip* in the other diagram, will therefore 
cancel each other in summations over all possible diagrams, as they contribute with equal 
amplitudes and opposite signs. As illustrated in Fig. IU this is the case for each diagram 
containing a block with ladder-type input, which has a counterpart in which this block is 
incorporated in the opposite manner. Such diagrams do therefore not need to be consid- 
ered for the calculation of mean densities or currents of the propagating condensate. Blocks 
with crossed-type input, on the other hand, can, in general, not be paired with canceling 
counterparts, which is shown in the right column of Fig. HI 

Let us finally point out that the validity of the present diagrammatic theory is still limited 
to weak nonlinearity strengths, despite the above summations to infinite order in g. This is 
illustrated in Fig. [5] which shows an example for a diagram of second order in g that is not 
accounted for in our diagrammatic theory. It represents diffraction of the matter wave by 
short-ranged spatial fluctuations of the nonlinear term 5f^|\l/(r,t)p in the Gross-Pitaevskii 
equation ([1]). As it requires the presence of two nonlinearity events within a region of the 
order of one wavelength, its contribution is strongly suppressed in the semiclassical regime 
as compared to other diagrams of second order in g in which the nonlinearity blocks are 
spatially uncorrelated. We do expect, however, that diagrams of the type shown in Fig. [5] 
will become relevant for large nonlinearity strengths, possibly in the regime in which the 
scattering process destabilizes and develops turbulent-like flow. 

4. Loop corrections 

In the previous section, we developed a semiclassical description of weak localization 
in the presence of a weak atom-atom interaction restricting ourselves to the diagonal ap- 
proximation. This theory will fail to describe the occurring phenomena quantitatively, as it 

31 



encounter region 




Figure 6: Sketch of a trajectory 7 which experiences a self encounter and its corresponding partner trajectory 
7'. In configuration space, which is depicted here, one trajectory crosses itself with a small crossing angle 
e whereas the other one avoids the crossing. Also shown is the position of a possible Poincare surface of 
section V used to determine the action difference. We note here that this sketch is widely overestimating 
the crossing angle e and underestimating the trajectory lengths before and after the encounter region. 



violates current conservation both in the absence and in the presence of the nonhnearity. The 
reason for this failure lies in the use of the diagonal approximation, i.e. we only used iden- 
tical or time-reversed trajectories when our methods demanded correlated trajectory pairs. 



However, as originally shown in Refs. |54j . |55| for the spectral form factor and in Ref. [30| for 
the Landauer conductance in the transport of electrons through two-dimensional uniformly 
hyperbolic systems with time-reversal symmetry, there is, besides the above mentioned one, 
a second type of correlated trajectory pairs giving significant contributions to the reflection 
and transmission probabilities, the so-called "loops" or "Sieber-Richter pairs". These are 
pairs of trajectories with nearly identical initial and final conditions; as illustrated in Fig. El 
one of the two trajectories undergoes a self-crossing with a small crossing angle e whereas 
the other one avoids that crossing. More generally speaking, as originally worked out in 



Refs. [56|, l57| for the spectral form factor, these trajectories exhibit an encounter in phase 
space with their time-reversed counterparts, which allows for the existence of a partner tra- 
jectory that switches from the original trajectory to the time-reversed counterpart. We shall 
adopt this more general phase space picture to derive the corrections to weak localization 
in the linear and in the nonlinear case. 

4-1. Loop corrections in the linear case 

Our calculation of the contributions to weak localization in the case (7 = mainly follows 



Refs. [58|, |53|. We shall restrict ourselves here and in the following to the case of at most 
one Sieber-Richter pair per diagram, as the presence of more such pairs would only result 
in higher-order contributions in td/th = l/{Nc + Nc). For the sake of definiteness, we 
shall consider the problem of determining the transmission probability tm = ja/j^ that is 
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associated with the scattering process from the incident channel i in the left lead to the final 
channel n in the right lead. Our purpose is therefore to calculate 
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As this quantity involves a product of two Green functions, we are concerned with sums 
over pairs of classical trajectories 7, 7' here. In the context of the diagonal approximation, 
we already evaluated in Section [3] the most dominant contribution to this transmission 
probability, for which 7' is identical to 7. 

The next group of systematically correlated trajectories consists of pairs 7, 7' that ex- 



hibit, as sketched in Fig. El a self-encounter in phase space [5J, l30|, |53| . Their action difference 



can be determined by defining a Poincare surface of section V within the encounter region, 
which is oriented perpendicular to 7 on the first passage of this trajectory through it, i.e., 
which is pierced by the first stretch of 7 at its origin. Linearizing the classical dynamics 
in the vicinity of this trajectory, we can define two basis vectors e^ and e„ within the two- 
dimensional surface of section V that are respectively oriented along the stable and unstable 
manifold of 7. The action difference between 7 and 7' is then evaluated as 
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su 



(130) 



where s and u denote the coordinates with respect to the basis vectors e^ and e^, respectively, 
at which the trajectory 7 pierces through V for the second time. Obviously, AS'-^^^/ can be 
sufficiently small, i.e. of the order of ^, if, as depicted in Fig. |6l one of the two trajectories 
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exhibits a self-crossing in configuration space with a very small crossing angle e 
The partner trajectory, whose existence and uniqueness is granted by the chaoticity of the 
classical dynamics, will then avoid that self-crossing and follow the loop in between the two 
piercings through V in the opposite direction. 

In order to evaluate the contributions of such Sieber-Richter pairs to Eq. fll29p . we 
need to determine the probability of a trajectory 7 to exhibit a near-encounter in phase 
space. Due to ergodicity, the probability density for the trajectory 7 to pierce again through 
the Poincare surface of section in the opposite direction at given coordinates s and u and 
after a given propagation time t2 after the first piercing is given by the Liouville measure 

— i/o(P) q)]/S(/i) with (p, q) the coordinates of the second piercing in the full phase 
space and S(/x)= J dPq' J (Pp'6[^ — Hq{p', q')] the phase-space volume of the energy shell. If 
we want to calculate the probability density for a trajectory 7 with a given total propagation 
time T to have a partner trajectory 7' with a given action difference AS^y^y = AS, we are 
tempted to integrate this Liouville measure over all "intermediate" propagation times t2 
between the first and the second intersection through the Poincare surface of section V, over 
all "initial" propagation times ti from the incident channel to the first intersection through V, 
over all "final" propagation times ^3 from the second intersection through V to the outgoing 
channel, as well as over all possible phase-space coordinates s, u that 7 exhibits within V 
at its second piercing, with the requirements that su = AS and ti + ^2 + i^s = T. This naive 
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integration would, however, lead to multiple countings of such trajectory pairs. Indeed, the 
placement of the Poincare surface of section V is not unique, but can be shifted along the 
first stretch of the trajectory 7. This generally will lead to different coordinates s, u of the 
second stretch of 7 when passing through "P, but the product su of these coordinates will not 
change, provided the second piercing point of 7 is also in a sufficiently close neighborhood 
of the origin of V such that the linearization of the classical dynamics around 7 is still valid 



(see also the calculations in Appendix E ) . 

The contribution of an individual Sieber-Richter pair with an action difference AS" would, 
when performing the above-mentioned integration, therefore effectively be overweighted by 
a factor tenc = ^enc(AS') that corresponds to the typical "duration" of the encounter, i.e., the 
typical propagation time within which one of the trajectories "sees" the other one within a 
distance that is within the linearization region of its transverse dynamics. Defining by A the 
Lyapunov exponent of the ergodic system 59|], and introducing c as the maximal distance 
along the stable and unstable manifolds e^, e„ for the linearization of the transverse dynamics 
within the Poincare surface of section V to be valid (i.e., we require that — c < s,m < c; the 
precise value of c, which is related to the Ehrenfest time of the system as pointed out in 



Appendix E.l will not be of relevance in the end), we can define [58|, |53 
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This duration tenc{su) refiects the fact that some minimal time is needed for the two nearby 
trajectory stretches to part from each other, in order to form the loop on one end and to 
exit toward different leads on the other end of the encounter region. 
In view of these considerations, we define (see also Ref. (53|) 
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as the probability density for a trajectory 7 with the total propagation time T to have a 
partner trajectory 7' with an action difference AS* and a loop duration ^2- This loop duration 
t2 as well as the initial and final propagation times ti and ^3 that appear in the integrations in 
Eq. fll32p are, as illustrated in Fig. [HI defined not with respect to the particular placement of 
the Poincare surface of section, but with respect to the location of the encounter region along 

the trajectory. Using y4^~y4^/ and /i^~7i^/ for the trajectory pair 7, 7', the loop contributions 

to {\ipf^\'^) are calculated as 
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where the dephasing factor exp(— )f:2/TB) originates from ^p^ [r7(ti + ^enc + t^)-, "s^-^itx + tenc), Z^], 
the flux integral along the loop according to Eq. (l52l) . We neglect here the contribution of 
the flux inside the encounter region, which will be discussed in the next subsection. 

When applying the sum rule (IC.lSp . we have to use a modified survival probability 
exp [— (T — tenc)/T"(i] in Eq. (lC.6p . Indeed, if the first stretch of the encounter lies within the 
billiard, the second one does so as well, thus the trajectory does not risk to escape during its 
second passing through the encounter region. The relevant time for the survival probability 
has therefore to be reduced by the duration tenc of this second stretch. We then obtain 
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td/tb ' 

as shown in Appendix E.l In a similar way as for nonlinearity blocks (see the discussion 
in Section l33|) . the encounter region cuts the diagram into three "links" contributing either 
Tn or td/{1 + td/tb). 

The same derivation can be applied in order to calculate the loop contributions to the 
reflection probability into channel n, leading to exactly the same result as in Eq. (11341) . We 
therefore obtain 
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as loop contributions to the reflection and transmission probabilities r^^^ and t^^} in the 
absence of interaction. These corrections do indeed restore current conservation in leading 
semiclassical order. Combining all ladder, crossed, and loop contributions that are evaluated 
in Eqs. (jl^ . f lFTI) . and fll35p . respectively, we obtain with Eq. 
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for the probability of retro-reflection into the incident channel n = i, as well as 
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for the probabilities of reflection into a different channel n ^ i and of transmission into 
channel h. This yields the total reflection and transmission probabilities 
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in the linear case g = 0, which obviously satisfy R^^'^ + T^^^ = 1. 
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Figure 7: Diagrammatic representation of the nonlinear loop contributions that arise in first order in the 
nonlinearity strength g. These diagrams involve encounters between different trajectories rather than a 
self-encounter of one trajectory. The upper and lower rows show, respectively, the two diagrams (a) and 

(b) in which the nonlinearity block is located outside the encounter region, as well as the three possibilities 

(c) , (d), and (e) for the nonlinearity event to enter the encounter region. The diagram (c), in which the 
nonlinearity block moves along a stretch through the whole encounter region, corresponds to the transition 
from diagram (a) to diagram (b). Diagrams (d) and (e), on the other hand, are obtained from diagrams (a) 
and (b), respectively, by pushing, in these latter two diagrams, the nonlinearity block along the trajectory 
that starts at the block into the encounter region. 

4-2. Contributions of first order in the nonlinearity 

In the case of nonvanishing interaction between the atoms, the determination of the 
loop contributions to reflection and transmission probabilities becomes more involved due to 
richer possibilities for associating correlated trajectories that exhibit small action differences. 
Loop contributions arise not only from self-encounters of single trajectories, but also from 
encounters of different trajectories in phase space. This is illustrated in Fig. [7] which shows 
the nonlinear diagrams that contribute to loop corrections of the reflection and transmission 
probabilities in linear order in g. As it is quite instructive, we begin our analysis of loop 
corrections in the nonlinear case with the calculation of the contributions of these diagrams. 
We shall first focus on diagrams (a) and (b) of Fig. [7] in which the nonlinearity event can 
only move along parts of trajectories that are outside the encounter region. 

As a starting point, we have to define the probability density w for having a pair of 
trajectories 7 and 7' that exhibit a near-encounter in phase space. This near-encounter 
results in the existence of an additional pair of partner trajectories7, 7'. In configuration 
space, the trajectories cross each other under a small angle in one of these two pairs, (7, 7') 
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(a) 



(b) 



(c) 



Figure 8: Three different constellations for the stretches involved in an encounter region. In the case of 
constellation (a), which arises also in the linear case discussed in Section I4TT1 the encounter region itself does 
not provide a significant contribution to the dephasing in the presence of an external magnetic field, since 
the flux integrals in the phases associated with the individual stretches cancel each other. This is different 
for the constellations (b) and (c) in which twice (b) and four times (c) the flux integral from a stretch of 
duration tcnc remains, yielding the phase factors exp(— tonc/'''B) for constellation (b) and exp(— 2ionc/TB) for 
constellation (c). 

or (7, 7'), whereas the other pair avoids this crossing. 

The probabihty density w is specified for a given action difference AS between the trajec- 
tories 7, 7' and the pair of partner trajectories 7,7', as well as for given partial propagation 
times t, t' of the trajectories 7 and 7', respectively, before (in the case of trajectory 7) or 
after (in the case of trajectory 7') the encounter region, which may become relevant for 
the evaluation of magnetic dephasing. It furthermore depends parametrically on the total 
propagation times T and T' of the two trajectories 7, 7' as well as on the orientations of the 
individual trajectory stretches within the encounter region, as these orientations might give 
rise to additional contributions to the magnetic dephasing. Figure [H] displays three different 
possibilities for orienting the trajectory stretches within the encounter region. While the en- 
counter region in constellation (a) does not contribute to the dephasing in the presence of a 
magnetic field, the constellations (b) and (c) contribute with phase factors exp(— tenc/'7"B) and 

exp(— 2tenc/TB), respectively, as there are one (b) and two (c) trajectory stretches ► 

of the linear Green function that are not balanced by the complex conjugate counterparts 



Denoting by n G {0,1,2} the number of imbalanced pairs of stretches within the en- 
counter region, we define 



as the density of trajectory pairs 7, 7' that come close to each other in phase space and 
thus have partner trajectories with a combined action difference AS. In this expression, the 
integration variables t and t' correspond to the propagation times of the final and initial 
parts of the trajectories 7 and 7', respectively, after leaving [t) and before entering [t') the 
encounter region. We have n = = Ua for diagram (a) and n = 1 = nb for diagram (b) in 
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Fig. El 
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Calculating now the contributions of the diagrams shown in the upper row of Fig. [7] 
(which are multiplied by a combinatorial factor 2 as there are two possibilities to construct 
these diagrams), we obtain 
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where we applied the sum rules (1C.14| IC.ISP to convert the sums over classical trajectories 
7, 7' into integrations over trajectory durations Ty, T^/, taking into account that we have 
to use a reduced effective time + — tcnc for the classical survival probability. Gauge 
invariance of the result is ensured by the fact that the encounter region closes the overall flux 
integral. The integration over s and u is calculated in Appendix E.l and yields —1/{tdTh) 
for n = na = as well as — (1 + td/tb) I{tdTh) for = = 1. We therefore obtain 
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for the cases where the nonlinearity block is located outside the encounter region. 

As shown in the lower row of Fig. [TJ there are two qualitatively different possibilities 
for the nonlinearity event to enter the encounter region. In the first scenario, depicted in 
Fig. [Zl^c), the nonlinearity event moves along a trajectory that provides a stretch within the 
encounter region. This scenario corresponds to the transition from diagram (a) to diagram 
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(b). Its contribution is calculated as 
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[see Fig. ^c) for the signification of ti, ^2, ^i, ^2]; where the integrations over s, u and over 
the propagation time t' within the encounter region (whose gauge field dependence is taken 
into account in the integration) yield, as shown in Appendix E.2 1/th- 

In the other scenario, depicted in Fig. UlA) and (e), a trajectory pair leaving the encounter 
and ending at a nonlinearity event becomes arbitrarily small until finally the nonlinearity 
event enters the encounter region but does not traverse it. This case requires a modification 
of the probability density of suitable partner trajectories as some stretches do not leave the 
encounter region any more but terminate at a certain point within it. Following Refs. 60|, 



6ll . |62| , we define a reduced encounter region with the duration 
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where t' G [0, (1/A) ln(c/|s|)] is the time interval between the nonlinearity event and the 
Poincare surface of section V that is optimally chosen to be located in the center of the 
encounter region. As a consequence, we have to extend the integration over s and u, asso- 
ciated with the possible action differences su, by an integration over all possible time spans 
t' defining the location of the nonlinearity event with respect to V, which substitutes one of 
the integrations over time in Eq. f ll40p . This yields the modified probability density 
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Using this density for the calculation of diagram (d) in Fig. [71 we obtain 
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where we evaluate the integration over s, u and ti in Appendix E.2 yielding 1/th- The 
calculation of the contribution of diagram (e) in Fig. [7] proceeds in perfect analogy with the 
one presented for diagram (d) and yields the same result 
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The overall nonlinear loop contribution to (^i°''*(5'?/'i^'')) originating from the diagrams 
shown in Fig. [7] therefore sums up to 
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As this expression is purely imaginary, no modifications of transmission and reflection prob- 
abilities are expected in linear order in the interaction strength g, which is in perfect accor- 
dance with the discussion in Section 13.41 [see Eq. f ll02p . 

4-.3. Contributions of arbitrary order in the nonlinearity 

As we can see from the above calculations, the presence of an encounter region perfectly 
flts to the picture of diagrams consisting of separated parts, which was developed in Section 
13.51 since we can perform the integrations corresponding to an encounter region indepen- 
dently from the remaining integrations over link durations. Encounter regions can, in the 
spirit of Section 13.51 be interpreted as extended "blocks" which are connected via four links 
to other (nonlinearity or encounter) blocks as well as to the leads of the system. Care must 
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be taken, though, if a nonhnearity block enters the encounter region, as is the case in the 
diagrams depicted in the lower row of Fig. [71 Under consideration of the calculations per- 



formed in Appendix E , we can extend our diagrammatic rules listed in Section 13.51 by the 



following ones: 

(7) each encounter region containing no nonlinearity event contributes a factor — (1 + 
nTB / Td) / {jdTh) where n = 0, 1, 2 counts the number of trajectory pairs with im- 
balanced stretches within the encounter region (see Fig. [H]); 

(8) each encounter region including a nonlinearity event contributes a factor 1/th- 



As was already argued in Section 13. 5[ diagrams containing a nonlinearity block with ladder- 
type input (as the ones shown in the left and central column of Fig. H]) will not contribute to 
the reflection and transmission probabilities as they are canceled by counterparts in which 
this nonlinearity block is attached to the complex conjugate trajectory stretch. We can 
therefore restrict ourselves to "crossed type" nonlinearity blocks shown in the right column 
in Fig. H 

Figure |9] shows the complete set of nonlinear diagrams contributing to reflection and 
transmission probabilities in arbitrary order in the nonlinearity strength g. These diagrams 
consist of the same chains of nonlinearity blocks as the ones shown in Fig. [31 with the 
main difference that these chains cannot be directly attached to a lead as we are calculating 
the scattering amplitude to an arbitrary (i.e., not necessarily the incident) channel in the 
reflected or transmitted lead; instead, they have to be connected to ladder links via an 
encounter region. In perfect analogy with Fig. [3 the nonlinearity blocks at the right ends 
of the chains in the diagrams (a), (b), (d), and (e) can be moved all the way through the 
entire encounter region, thereby giving rise to a transition from diagram (a) to diagram (b) 
as well as from diagram (d) to diagram (e). The other blocks located at the ends of the 
chains in the diagrams (a) - (g) can, in analogy with the diagrams (d) and (e) of Fig. [71 be 
pushed into the encounter region by moving them along the trajectories that start or end 
at those blocks. 



Using the diagrammatic rules (1) - (8) stated in Section [3751 and above, the contributions 
of these relevant diagrams as well as their corrections due to nonlinearity events entering 
the encounter region are straightforwardly evaluated. Defining by 

90 - (/r,4»') = hr,'JL^_\ " (150) 



41 




Figure 9: Complete set of diagrams of k-ih order in g containing one encounter region and k nonlinearity 
blocks of the "crossed" type shown in the right column of Fig. U] which are located outside the encounter 
region. The diagrams (a) - (f) are contributions coming from (V'ff^*('5V-'i'^^)) ^nd (('^V'^f respectively, 
where (c) and (f) only exist for k > 2. Contrary to the other diagrams, (g) is contained in ((<54'''^*)(^V'i'''^)) 
with ki + k2 = k and fci,fc2 > 1. As is shown in Eq. (I158|) . the contributions of the diagrams (b), (c), (e), 
and (f) exactly cancel each other. 
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the relevant scale for the nonlinearity strength appearing in Eq. (11241) . we obtain for the 
diagrams (a) - (g) 
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where dif)^^^'^^'^^ denotes the /c-th order contribution to ipn, whose diagrammatic representa- 
tions contain ki nonlinearity blocks of the type I and k2 complex conjugate nonlinearity 
blocks H . The corrections to these contributions due to nonlinearity events entering the 
encounter region are calculated as 
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By multiplying these individual contributions with the corresponding powers of g and 
then summing over all orders k, ki and ^2 in analogy with Eq. (11281) . we finally obtain 
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for the summed contributions of the diagrams (a) and (d) as well as for the contributions of 
the diagrams (g), and 

/I / |2\ |2\ _ o mf frpV {g/gof 

/loop, (b) + (e) - ;i„„p, (,) + (f) - y^^j ^ ^ ^^^^^^2 U^SJ 

for the summed contributions of the diagrams (b) and (e) as well as of (c) and (f), which 
implies that the contributions of these latter four diagrams exactly cancel each other. The 
summation of the associated corrections due to nonlinearity blocks entering the encounter 
region yields the contributions 

M'''"! /loop, (a) +(d), correct \ I '''"I /loop, (c) + (f), correct p\{fi) \Th J I + {q / QqY ' 
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which exactly cancel each other as well. In effect, therefore, only the diagrams (a), (d), and 
(g) provide nonvanishing contributions to the reflection or transmission probabilities, which 
are summed up as 

<IV^^I'>ioop,{a)+(d)+(g) = i + td/tb {£) 1 + {g/goY ■ ^^^^^ 

Together with the result obtained in the linear case, the overall loop contribution to 
(IV'nP) reads 
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which yields as the correction to the transmission and reflection probabilities 
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As in the linear case [see Eq. f ll35p ]. this correction restores current conservation in the 
presence of the nonlinearity. With Eq. fll24p we obtain 
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for the probability of retro-reflection into the incident channel n = i, as well as 



N, + N, {N, + N,)Hl + rn/TBf + {gfrof 

for the probability r„j of reflection into a different channel n ^ i oi the incident lead and 
for the probability tni of transmission into channel h. This yields the total reflection and 
transmission probabilities 
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which obviously satisfy R + T = 1. 

5. Comparison with numerical results 

Figures [10] and [11] display (in their right columns) the semiclassical prediction fll64p for 
the probability of retro-reflection into the incident channel as evaluated for the billiards 
a and b, respectively, that are shown in Fig. [T] The sizes of the two billiards are chosen 
such that both the incident and the transmitted leads exhibit five open channels, i.e. Nc = 
Nc = 5, at the energy that corresponds to the chemical potential fi of the incident beam. We 
specifically have the areas ~ 3.41 x 10^ h^/ (m/xo) for bilhard a and ~ 3.29 x 10^ h^/ (m/io) 
for billiard b, where fiQ defines the characteristic energy scale for the chemical potential 
of the atomic beam (i.e. we choose fi = Hq for the evaluation of the semiclassical retro- 
reflection probability). The incident current is chosen as f = l.Ofio/H. As described in 



Appendix D , the dwell time To and the characteristic scale Bq of the effective magnetic field 
were classically determined from the numerically computed length and area distributions 
within the two billiards, respectively; we obtained fro ~ 267 and Bq ~ 1.55 x 10~^ mfiQ/h ~ 
0.844 X 2TTh/n for billiard a as well as jV^ ~ 241 and Bq ~ 4.21 x 10-^m^o/^ - 0.221 x 
2TTh/Q for billiard b. 

In the linear case g = 0, a. Lorentzian peak is obtained for the retro- reflection probability 
as a function of the effective magnetic field, on top of an incoherent background at ra ~ 
l/{Nc + Nc) = 0.1. This is the characteristic signature of weak localization. As is evident 
from Eq. fll64p , the presence of a finite nonlinearity g gives rise to a reduction of the coherent 
enhancement of the backscattering probability, which ultimately approaches the incoherent 
background l/(A'^c + A'^c) for g — t- oo. This reduction, however, is more effective at the center 
of the backscattering peak than in its wings, such that for intermediate values of g a local 
minimum may be encountered in the reflection probability around B = 0. 

This prediction is indeed confirmed by numerical computations of the quasi-stationary 
transport process within the two billiards under consideration. As explicit numerical prop- 
agations of the time-dependent Gross-Pitaevskii equation are rather time-consuming, 
we use, in practice, a different approach in order to compute the scattering states of the 
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Figure 10: Nonlinearity-induced destruction of weak localization for the billiard a displayed in the left panel 
of Fig. [TJ Plotted are the numerically computed backscattering probabilities (upper left panel) and their 
semiclassical prediction according to Eq. (|164l) (upper right panel) as a function of the effective magnetic 
field for various values of the nonlinearity g. We use the magnetic field scale Bq ~ 1.55 x 10~^m/zo/?i and 
the average population j'td ~ 267, which were inferred from an analysis of the classical dynamics within the 
billiard. The lower panels display the differences of the backscattering probabilities for finite g with respect 
to the backscattering probabilities of the linear system. Good agreement is found between the numerical 
data (lower left panel) and the semiclassical prediction (right panel). 
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Figure 11: Same as Fig. [10] for the billiard b displayed in the right panel of Fig. [T] We use the magnetic 
field scale Bg ~ 4.21 x 10~'^ mij^o/h and the average population jWu ~ 241, which were inferred from an 
analysis of the classical dynamics within the billiard. 
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Figure 12: Total transmission for the billiard a. Plotted are the numerically computed transmission proba- 
bilities (upper left panel) and their semiclassical prediction according to Eq. (jl66|) (upper right panel) as a 
function of the effective magnetic field for various values of the nonlinearity g, using the same parameters as 
for Fig. 1101 The lower panels display the differences of the transmission probabilities for finite g with respect 
to the transmission probabilities of the linear system. Good agreement is found between the numerical data 
(lower left panel) and the semiclassical prediction (right panel) . 



48 



1 ^ \ ^ \ ^ \ ^ f 




Figure 13: Same as Fig. [T^] for the billiard b. 
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system. As explained in Appendix A this approach uses a Newton search algorithm in 
order to construct self-consistent scattering states of the nonlinear system. Among all those 
scattering states that are identified at given chemical potential fj, and given incident current 
(there is only one such scattering state in the linear system, but several of them generally 



exist in the presence of the nonlinearity, see, e.g., Ref. [63[), we select the one that would 
be first encountered in the presence of an adiabatic increase of the nonlinearity strength g. 
More technical details concerning this approach will be provided in a subsequent publication 
0. 

The left columns of Figs. [10] and [11] display the results that are obtained from these 
numerical computations. To obtain good statistics, we did not only perform an energy 
average of the reflection probability, within the energy interval 0.93 /xq < /i < 1.18 /io for 
which there are exactly 5 open channels within each lead, but also averaged over different 
positions of the semicircular and circular obstacles in the case of billiard a and b, respectively. 
This additional configurational average is necessary as the above energy interval contains 
only a limited number of resonances within the billiard. Moreover, we averaged over different 
choices of the incident channel i, even though only the choice i = 1 appears realistic from 
the experimental point of view. The error bars attached to the data points consequently 
indicate the size of the statistical standard deviation that results from these averages. 

As shown in the upper panels of Figs. [10] and [Til the relative height of the peak with 
respect to the incoherent background significantly deviates from the universal semiclassical 
prediction, even in the linear case g = 0. This discrepancy may, on the one hand, be 
attributed to a limited applicability of the semiclassical framework in our context. Indeed, 
as is seen in Fig. [T] the wavelength of the matter-wave beam is not sufficiently small to rule 
out the influence of possible diffraction effects at the rounded corners of the billiard. On 
the other hand, non-universal scattering phenomena that explicitly depend on the shape 
of the billiard under consideration may play a role. Specifically, among the backrefiected 
trajectories that start and end in a given channel, there is possibly a significant fraction of 
self-retracing trajectories which are identical to their time-reversed counterpart. As those 
self- retracing trajectories do not contribute to the crossed part of the coherent backscattering 
probability, their semiclassical contribution then should be subtracted from the sum-rule 
based expression dSSj) of the crossed density. Indeed, the presence of a prominent circular 
obstacle within billiard b should allow for a number of rather short self-retracing trajectories 
with a relatively small Lyapunov exponent (and therefore with a relatively large weight in 
the semiclassical Green function), namely those trajectories that directly head toward the 
obstacle, undergo a self-retracing reflection there, and subsequently exit the billiard in the 
incident channel. Similarly relevant self-retracing trajectories bouncing off the semicircular 
obstacles should exist in billiard a. 

In view of the diagrammatic theory developed in Sec. [3] we note that such self-retracing 
trajectories do not affect the nonlinearity- induced corrections c^f — c^^^ to the coherent 



backscattering probability. Indeed, as is evident e.g. from Eq. fllOSp . those corrections are 
distinctly different from ladder contributions and will therefore not be overcounted if they 
happen to involve self-retracing trajectories. We consequently find, as shown in the lower 
panels of Figs. [10] and [HI rather good agreement between the numerical data and the 
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Figure 14: Change of height of the weak locaUzation peak (or dip) at B = as compared to the hnear 
case. Plotted are, in the left panel, \rii[g = 0) — rii(^g^{Nc + Nc)"^ /{Nc + Nc + '^) and, in the right panel, 
\T{g) — T{g — 0)] {Nc + /N^ as a function of gfro / {Nc + Nc) for the two biUiards under consideration. 



A comparison with the semiclassical predictions (|164l) and (|167p (solid lines), which are given by 1/(1 - 
with X = gfTujiJ^c + -^c) in both cases, shows good agreement. 







semiclassical prediction if we specifically compare those corrections, i.e. the reduction of the 
weak localization peak with respect to the linear case g = This is furthermore confirmed 
in the left panel of Fig. [H] which shows the reduction of the height of the weak localization 
peak at S = as a function of the nonlinearity strength g. Renormalizing the horizontal and 
vertical axes in terms of the scales that are suggested by the analytical prediction f ll64p . we 
find rather good agreement with this universal prediction for both billiards. This underlines 
the validity of the approach developed in Sections E] and IH 

Finally, in order to demonstrate the relevance of the loop contributions, we show in 
Figs. [121 [121 and in the right panel of Fig. [TDthe comparison of the numerical results for the 
transmission with our analytical prediction obtained from Eq. fll67p . Once again, excellent 
agreement is found after removing non-universal effects. Remarkably, Figs. [T2] and [T5] display 
asymmetries in the transmission as a function of the magnetic field, i.e. we do not necessarily 



have T{—B) = T{B). This finding seems to constitute a violation of Onsager's relations [65 
which state that the total mesoscopic transmission T = J2i I'^-niP; which represents the 
incoherent sum over the individual transmissions that result from all available choices of the 
incident channel i (and which is implicitly computed in Figs. [T^ and [T^ due to the averaging 
over the incident channel), be symmetric in the magnetic field B for any given scattering 
geometry at any given chemical potential. It should be noted, however, that Onsager's 
relations are based on the unitarity of the scattering matrix S and its symmetry property 
S{—B) = ^"^(i?) [66] and thereby implicitly rely on the linearity of the scattering process. 
Indeed, computing the total mesoscopic transmission T across billiard b in the absence of 
interaction for a specific choice of the chemical potential, we obtain perfect symmetry of T 
in B as shown in Fig. [151 This symmetry is broken at finite values of the nonlinearity g. 
Similar findings have been reported in electronic transport through mesoscopic structures 
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Figure 15: Total mesoscopic transmission T = ^ X^iX^fi l^niP a function of the magnetic field B for 
billiard b at a fixed chemical potential n = 1.08306 fiQ and at different strengths of the nonlinearity g. In the 
absence of interaction (black curve), T is symmetric in B, which is a consequence of Onsager's relations (sHj . 
This symmetry in B is, however, broken in the presence of the nonlinearity (red, blue, and green curves). 
A sufficiently large energetic and configurational average would restore the symmetry in B as is visible in 

Fig.m 




in the presence of strong bias voltages 67|, |68 



6. Conclusion 

In summary, we studied, both analytically and numerically, weak localization of guided 
matter waves that originate from interacting Bose-Einstein condensates and propagate through 
chaotic billiard geometries. Our analytical approach is based on a nonlinear diagrammatic 



perturbation theory [32|, |33|, |3J] that originates from the Gross-Pitaevskii equation, which 
is combined with a semiclassical expansion of the linear (single-particle) Green function 
within the billiard. Summing all terms of this diagrammatic perturbation theory and utiliz- 
ing standard semiclassical sum rules in ergodic billiards, we obtain analytical expressions for 
the retro-reflection probability [Eq. fll64p ] as well as for the total reflection and transmission 
[Eqs. fll66p and fll67p ] in dependence of the effective interaction strength and of the strength 
of an artificial gauge field that breaks time-reversal invariance and simulates the effect of 
a magnetic field for charged particles. These expressions also involve the analysis of loop 
corrections in leading order [30| which restore current conservation. 

Globally, we find that the peak of weak localization decreases with increasing nonlinearity 
strength g and eventually disappears beyond a characteristic scale of g given by the inverse 
average population within the billiard. This suggests that the presence of the nonlinearity 
introduces an additional dephasing mechanism that affects the constructive interference 
between backscattered trajectories and their time-reversed counterparts. The decrease of 
the peak height with g is found to be stronger at the center of the peak (i.e. for vanishing 
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gauge field -8 = 0) than in its wings, which eventually gives rise to a tiny local dip in the 
backscattering probability around B = 0. While this dip, as it is predicted by the general 
semiclassical theory, is presumably too small to be of experimental relevance, it can be more 
pronounced in individual billiard geometries, as the ones specifically studied in this work, in 
which the backscattering probability develops a global minimum, instead of a maximum, at 
B = 0. We thereby encounter a signature of weak antilocalization in those billiards, which is 
of genuinely different nature than antilocalization in electronic transport processes involving 



spin-orbit interaction [28 



Comparisons of the numerically computed absolute and relative heights of the weak local- 
ization peaks with the semiclassical prediction seem to suggest that this weak antilocalization- 
type phenomenon is caused by the occurrence of self-retracing trajectories in the scattering 
system. Indeed, the presence of such self-retracing trajectories reduces the probability for 
coherent backscattering as compared to the universal semiclassical prediction in the linear 
case, as the application of the standard sum rule would give rise to an overcounting of 
interference contributions between backscattered trajectories and their time-reversed coun- 
terparts. It does, however, not affect the nonlinearity-induced corrections to this coherent 
backscattering probability. Consequently, the peak of weak localization can turn into a finite 
dip in billiard geometries that exhibit prominent self- retracing trajectories of short length 
and therefore of large weight in the semiclassical backscattering amplitude. 

This observation also sheds new light on the inversion of the coherent backscattering 
peak that was found in the coherent propagation of Bose-Einstein condensates through 



two-dimensional disorder potentials [27[. As a matter of fact, such disorder potentials also 
exhibit self-retracing trajectories, which essentially arise from a retro-reflection at the first 
impurity that the incident matter wave encounters within the disorder region. Diagrammatic 
calculations within such disordered systems ^] do indeed suggest that short reflected paths 
are at the origin of the inversion of the coherent backscattering peak in disordered systems. 

In this study, we considered a number of idealizations concerning the setup for the matter- 
wave transport process. For the sake of analytical tractability of the problem, we particularly 
imposed hard-wall boundaries of the wave guides and the billiard and assumed a continuous 
monochromatic flow of atoms through this scattering region. We are convinced, however, 
that the phenomenology studied in this work should be sufficiently robust to manifest also in 
the case of harmonic waveguides and harmonic-like confinement geometries with chaotic (or 
mixed regular-chaotic) dynamics, which could possibly be realized by combinations of red- 
and blue-detuned laser beams that are perpendicularly focused onto the waveguide 0] , as well 
as in the case of atomic wave packet scattering processes which may be easier to realize than 
guided atom-laser beams. Weak localization and antilocalization of interacting Bose-Einstein 
condensates should therefore be observable with present-day cold-atom technologies. 
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Appendix A. Numerical computation of stationary scattering states 



In this appendix, we explain how we numerically compute stationary solutions ip = 
ip{r) of the time-dependent inhomogeneous Gross-Pitaevskii equation ([5]). Such stationary 
solutions satisfy the nonlinear equation 



^?(r)^|^(r)r^(r) 



Sir) = 



(A.l) 



with S{r) = SoXi{y)S{x — xl), which is equivalent to Eq. f l62|) . This equation is discretized 
on a two-dimensional lattice where only points inside the cavity and the leads are taken into 
account. The single-particle Hamiltonian H given by Eq. can be approximated using a 
finite-difference scheme 65| where we incorporate the vector potential using a Peierls phase 
(7ol |. We choose the lattice spacing A small enough that the approximation error [which 
scales as C(A^)] becomes negligible, which is the case for roughly 30 lattice points per 



wavelength. The interaction strength (yf(r) is, as explained in Ref. 71 



constant within the scattering region and adiabatically ramped off in the leads [72 



considered to be 
The 



effects of the infinite leads can then be incorporated through self-energies as in the recursive 
Green function method 66|, |73| , which allows one to restrict the numerical computation to 
a finite spatial region. 

The complex solution ip{r) of the nonlinear wave equation flA.ll) can be represented as a 
2A/'-dimensional real vector where J\f is the number of grid points. Defining 



F : 



^(r) I — ^ - ^^o] ^(r) - ^(r 



^2 
— I 

2m' 



:r)rV^(r)-5(r) 



(A.2) 



we search now for a solution of -F(V') = 0. This is done with Newton's method [74|. One 
selects a starting vector '?/'o(r) and constructs a sequence of vectors {4'k}kLi (here k is the 
iteration number) using the iteration ^pk+l = V'fc ~ {1^F)~^ F{iljf^). If the derivate T>F at 
the solution is not singular, this iteration is guaranteed to converge to a solution of the 
nonlinear equation (]A.1|) . provided the starting vector is chosen in a suitably close vicinity 
of this solution. 

The efficiency of this method strongly depends on the starting vector iPo{t^)- An obvious 
choice would be the solution of the linear wave equation (for = 0). This choice, however, 
works out only for very small nonlinearities. In the general case, one has to use a continuation 
method 7^,l75|. To this end, we consider g, i.e. the constant value of the interaction strength 



within the billiard 
function F : M^A^x 



76 



as an additional free parameter and reinterpret F = F[ip{r)] g] as a 
— )■ M^-^. Now F~^(0) is a one-dimensional manifold ff?] which can be 
conveniently parametrized by the arclength s through the parametric curve s t— )■ [g{s), ip{s)]. 
Starting from g = 0, the numerical algorithm follows this curve until the desired value g of 
the nonlinearity strength is reached, and returns the wavefunction tp{r) that is obtained at 
the end of this curve-tracking process [78 . 
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Figure A. 16: Transmission of the stationary scattering state versus the interaction strength g as obtained 
from the numerical curve-tracking method. Black lines indicate dynamically stable branches and red lines 
indicate dynamically unstable branches. Due to the nonlinearity in the wave equation, several scattering 
states with different transmissions may be encountered at a given value of g. At g = 0.08 for instance, 
there are three different scattering states, two of them being dynamically stable (black dots) and one being 
unstable (red dot). The light-blue data points show the time-averaged transmissions that are obtained 
from a numerical simulation of the time-dependent scattering process in the presence of an adiabatically 
slow ramping of the source amplitude. As is visible from the error bars which indicate the associated 
standard deviation of the transmission at given value of g, the atomic current through the billiard develops 
a pronounced time dependence in the absence of stable stationary scattering states (e.g. around g = 0.25). 
The calculation is done for billiard b at £? = —0.001 m/io/ft and fi = 0.935102 hq, with the condensate being 
injected in the transverse ground mode {i = 1) of the incident lead. 



Fig. lA. 16) shows, for a specific set of parameters, a projection of this curve onto the two- 
dimensional parameter space spanned by the nonhnearity strength g and the total trans- 
mission that is associated with the stationary scattering state ipir). As is characteristic for 
nonlinear wave equations, several stationary solutions are found for some values of g, e.g. 
a.t g = 0.08. Some of these solutions may be dynamically unstable and can therefore not 
be populated in the time-dependent propagation process. At sufficently large values of the 
nonlinearity, no dynamically stable scattering state is found any longer, which implies that 
the scattering process becomes permanently time-dependent and develops turbulent-like be- 



haviour 43 



To determine the dynamical stability of the numerically computed stationary scatter- 
ing state, we linearize the time-dependent Gross-Pitaevskii equation around the stationary 
solution 'ip{r) using the Bogoliubov ansatz 

ip{r,t) = [ipi^r) + u{r) exp {—iC,t) + v*{r) exp {zCt)]exp(-'-fit] (A.3) 
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Figure A. 17: This figure shows as a function of g the fraction of configurations of B, n and the incident mode 
i that exhibit only dynamically unstable solutions (red curve) or support more than one dynamically stable 
solution (blue curve). In the latter case, we select the one that is first encountered in the curve-tracking 
algorithm starting from g = 0. 



for the time-dependent scattering wavefunction ilj{r,t). This leads to the Bogohubov-de 
Gennes equation 

'H-f^ + 2g{r)^mr)\' ^^Wfi [^(r)]^ \ M\ ^ (+1 \ (u{r) 

9ij)t \rir)f H*-f^ + 2^(r)£|^(r)p; ^^^(r); ^ -l) {vir) 

(A.4) 

This generahzed eigenvalue problem is numerically solved using the implicit restarted 
Arnoldi method as realized in the software library ARPACK [80]. Special care must be taken 
in order to describe the infinite leads properly. This is done using the method of smooth 



exterior complex scaling [81| which exponentially damps the outgoing waves of the collective 
modes in the leads. As a consequence, the Bogoliubov eigenfrequencies ^ become complex. 
If one of them is found to exhibit a positive imaginary part, i.e. Im(,^) > 0, we can infer 



from Eq. flA.Sp that the scattering state ip{r) under consideration is dynamically unstable. 



Fig. IA.17I shows (red curve) the fraction of parameter configurations of the chemical po- 
tential fi, the magnetic field B, and the incident mode number i for which no dynamically 
stable scattering state is found. This fraction of unstable configurations is found to increase 
with the nonlinearity strength g, which imposes restrictions on the shape of the cavities and 
the maximum value of g one can use for numerical simulations. In particular, we find that 
the fraction of configurations that support no stable solution increases rather rapidly with g 
for small widths of the leads, i.e. for a very low number of open channels. This is attributed 
to the reduced spectral width that quasi-bound resonance states within the billiard exhibit 
in that case, which in turn leads to an enhanced interaction energy J drg{r)\tp{r)\'^ at such 
resonances. Rather wide leads with a large number of open channels, on the other hand, 
compromise the effect of weak localization and reduce the visibility of its signature in the 
reflection and transmission probabilities. For the billiard sizes and geometries under con- 
sideration, the fraction of conflgurations with unstable solutions remains below 0.02 until 
g = 0.08. 

When encountering a conflguration with only unstable solutions, we select the one that 
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exhibits the smallest Lyapunov exponent, the latter being defined by the largest imaginary 
part of the eigenvalues ^ of the Bogohubov-de Gennes equations (1A.4|) . This choice is sup- 
ported by time-dependent propagations of the inhomogeneous Gross-Pitaevskii equation ([5]), 
which directly simulate, as in Ref. 27|, the time- dependent scattering process in the presence 
of an adiabatic increase of the source amplitude. In such simulations, which were carried out 
for specific parameter configurations, we find, for not too large values of g, that the time- 
dependent current in the transmitted lead displays regular oscillations around the current 
of the stationary solution that is the least unstable one and exhibits the smallest Lyapunov 
exponent. The time- averaged transmission, which is the main experimental observable in 
such scattering processes, is then correctly reproduced by the transmission of the unstable 
stationary scattering state. For larger values of g, however, the time-dependent dynam- 
ics of the propagating wavefunction becomes chaotic, which means that the transmission 
associated with an unstable stationary scattering state loses its significance. 



Appendix B. The eikonal approximation 

In this appendix, we explicitly derive the semiclassical expression for the Green function 
in the presence of a weak perturbation of the Hamiltonian. Considering the Hamiltonian 
H = Hq + SH with Hq = Hq{p, r) the unperturbed part and 6H = SH{p, r) a perturbation 
that is slowly varying with r, and defining Gq = {E — Hq + ^0)~^, we can express the Green 
function of the system by means of a Dyson equation of the form 

oo 

G = {E-H + iO)-^ =Go + GoSHG = Go Y.i^HGo)'' . (B.l) 

Let us first evaluate the first-order term of this Born series, 

5G^'\r, r', E) = {r\Go5HGoW) = J rfV'Golr, r", E)5Hip, r")Go(r", r', E) , (B.2) 

in the semiclassical approximation. Using the semiclassical expression fllUI) for the Green 
function, 

Go(r,r',E) = ^A^(r,r',E)exp 

7 

we see that the momentum operator p in 6H acts, in leading semiclassical order (i.e., in 
lowest order in h), only on the action integrals S^{r" ,r' , E) in the exponents. This means 
that p can be replaced by the final momenta of the trajectories leading from r' to r". We 
therefore obtain 

SG^'\r, r', E) = I dh" ^7.(r, r", E)A,,ir", r', E)SHip", r") 

71.72 

X exp |^[S^,(r, r", E) + S^,(r", r', E)]-t^ (^^, + /x^jj (B.4) 
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(B.3) 



with p" = pi^_^{r" , r' , E) , where the indices 71 and 72 respectively represent the trajectories 
from r' to r" and from r" to r. 

Using now the fact that A^-^, A^^, and 6H {p",r") are slowly varying functions of r" on 
the length scale of the de Broglie wavelength of the atoms, we can apply the stationary 
phase approximation to evaluate the integral over r". The stationary phase condition yields 



^[S,,{r,r",E) + S,,{T",r',E)] = 0. 



(B.5) 



i.e., p^^ (r", r', £■) = p!^^(r, r", i?). This condition is satisfied if and only if the trajectory 72 



72 



is the direct continuation of 71. The double sum in Eq. (]B.4p can therefore be contracted 
to a single sum over trajectories 7 that are going from r' to r at energy E. Combining the 
prefactors that result from the spatial integration perpendicular to this trajectory as well 
as from the amplitudes A^-^ , A^^ , and transforming the spatial integration parallel to the 
trajectory into an integration along the propagation time, we finally obtain 



6G^'\r,r',E) = A,(r, r', E)exp 



6H[p^{t),q^{t)]dt. 

(B.6) 



Similarly, higher order terms in the Born series (IB.ip can be evaluated yielding 
6G^''\r,r',E) = {r\GoiSHGo)''\r') 
= ^A'yir,r',E)exp 



1 TT 

-S^{r,r',E) -i-fi^ 



X- 



6H[p^{t),q^{t)]dt 



(B.7) 



This finally yields the modified Green function 

G(r,r',E) = 5^A,(r,r',i?)exp 

7 

where 



^S^{r,r',E) -i^fi^ 



>^{r,r',E) 



{r,T',E) 



SH[p^{t),c^.^{t)]dt. 



(B.8) 



(B.9) 



represents an effective modification of the action integral 5*^ due to the presence of the 
perturbation. Eq. f IB.Sp reflects a general result in classical mechanics that, to leading order 
in the perturbation, the action difference between unperturbed and perturbed orbits is, for 
periodic orbits, given by Eq. (IB.Qp [82|. 

In the case of a perturbation due to a weak magnetic field, we have 



SH{p,r) = --A{r).p+^A\r) 
m 2m 



(B.IO) 
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Hence, we can write 0^ = —ip^ — ip^ with 



1 



m 



7 



<^^(r,r',E) = -/ p^(t)- A[q^(t)](it, (B.ll) 







V5^(r,r',i?) = -1-^^" K\q_^{t)\dt (B.12) 

corresponding, respectively, to a paramagnetic and a diamagnetic contribution to the effec- 
tive action integral. 

Also the presence of a weak nonlinearity within the scattering system can be accounted 
for in this framework, provided only ladder contributions are considered. Comparing Eq. (166|) 
in the cases (168|) and f l69|) with Eq. ( 1B.2|) . we see that we have to set 

m{r) = 2g^\i,^^{j)\l (B.13) 

where |'?/'*^°-*(r)|^ represents, according to Eq. fl72l) . the density at position r as evaluated 
within the diagonal approximation. We then obtain the effective modification of the action 
integral [defined by x-y(r, r', \x) in Sec. 13. 3j as 

-2 /-T. 



r,r',/x) = <7- / '|^(°)[q,(t)]|^rft. (B.14) 
^ Jo 



Appendix C. Sum rules 



In this appendix, we derive the generalized Hannay-Ozorio de Almeida sum rules |48|, |49 
that we need in order to evaluate energy averages of squares of the Green function in the 
diagonal approximation. To keep the derivation as general as possible, we introduce a new 
parametrization of the initial and final phase space points according to (p, r) = {C,i,C,2, C,3, ^4) 
and (p', r') = {r]i,f]2,V3yV4) where the sets (^1,^25^35^4) (^15 ''?2 5 '735 '74) contain the com- 
ponents {px,Py, X, y) and {p'^^Py, x' , y') of the final and initial phase space points, respectively, 
in some arbitrary order. We shall now be interested in the Green function from the coordi- 
nates {j][,rj'2) to the coordinates (^1,^2)- In the diagonal approximation, the energy average 
of the modulus square of this Green function reads 

GUi,i2\{v[,i2\Ef)^ = $^(lA[(656)5(r]i5r]^)5i^]|') (C.l) 



7 



^[(^156)5^5^/^)5^?] 



(C.2) 
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We furthermore need the corresponding expression for the crossed average which includes, 
in addition, a magnetic dephasing. This yields 



{G* U, V2), (^1, 6), E] G [(6,6), V2), E]) 
= 5^(|A,[(6,6),(^i,^2),^]l')exp 

d[ip',,p'y),ix',y'),T] 



-y 



det 



d[{^i,^2),{v'i,V'2),E] 



exp 



TV] 

^b) 



(C.3) 
(C.4) 



where tb ~ B^^ [see Eq. ( 153|) ] is the dephasing time. 

Applying standard rules for multidimensional integrations over (5-distributions, we can 
derive 



E 



det 



dM,p'),{x\y'),T] 



dU,,i2),{i,,v'2),E] 



j d^3 j j dv's j dT^'^j^ dT6[E-Ho{p',r')] 
x6[r- q(p', r', T)] 6 [p - p(p', r', T)] /(T) (C.5) 



for any /(T), where we define q(p', r', T) = {q^, qy){p', r', T) and p(p', r', T) = (px,Py)(p', r', T) 
as the position and momentum variables that result from the propagation of a classical trajec- 
tory over time T with the initial values p' = {p'^iP'y) and r' = (x', y'). Furthermore, assuming 
classical ergodicity, which is valid if the dynamics within the billiard is fully chaotic, we can 
state that each phase-space point (p, r) within the billiard has equal probability to be hit by 
a given trajectory after a given evolution time T, provided it lies within the shell of constant 
energy E. This probability, however, decreases exponentially with the evolution time, due 
to the possibility for escape from the billiard via the openings. We therefore obtain 



(5[r-q(p',r',T)]5[p-p(p',r',T)]) 



^[go(p^rO-go(p,r)] (_T_ 
J d'p J d^qS[Hoip',v') - Hoip,cDf'''^ V rn 



(C.6) 

where the "dwell time" td corresponds to the mean evolution time that a classical trajectory 
spends within the billiard before escaping to one of the waveguides. This altogether yields 



\G[{^,,^,),{VvV2),E]\' 



TD I d^3 J d^,S [E - Hoip, r)] / dr]', J di,5 [E - /7o(p', r')] 



Jd^pfd'q6[E-Hoip,q)] 



(C.7) 



and 



(G* H, V2), (^1, 6), E] a [(6, 6), iv'i, V2), E]), 



1 + td/tb 



G[i^u^2),iv'i,V2),E]\- 



(C.8) 

The phase space integrals appearing in Eq. f lC.7p can be straightforwardly computed. 
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We obtain 



Jd\6[E-Ho{p,q)] = n6(^E~^y 
j (fp6 [E - Ho{p, q)] = 27rmxn(q) , 
(fp j (fqS[E - Ho{p,q)] = 2-KmVt , 



(C.9) 
(C.IO) 
(C.ll) 



where Xn(q) represents the characteristic function of the scattering system and VL denotes 
the area of the bilhard. Furthermore, for the case of "mixed" initial or final conditions 
z = {xL,Py) specified within the incident lead, we calculate 



dy j dp J [E-i/o(p,q)] = ^ 



W poo 

dy J dpj ( E 



pI+pI 



mW 



(C.12) 



Here the longitudinal momentum p^ is restricted to positive (or negative) values correspond- 
ing to an initial (or final) condition with an incoming (or outgoing) trajectory. The width 
W of the waveguide is to be replaced by W in the case of a final condition within the 
transmitted lead. 

Putting these ingredients together and specifying the choice of the phase space variables 
(^15^2) and {rii,ri2) that appear as arguments in the Green function, we finally obtain 



|G(r,z',E)| 



G{z,r,E) 



G(z,z',E) 



TH\hy ^''^'"^271 ^2mE-p'^ 



Td fm\- 



W 



TH\hV ^''^'"^27r ^2mE-pl 
To /m\2 WW 1 



TH\h?J (27r)2 ^2mE - pf^ ^2mE - p 



(C.13) 
(C.14) 
(C.15) 



for z, z' being defined in the transmitted and incident lead, respectively, where th = mQ/h 
denotes the Heisenberg time of the billiard. The corresponding energy-averaged crossed 
densities are obtained by a multiplication with the prefactor (1 + to/tb)~^, as is seen from 
Eq. dUS}. 



Appendix D. Analysis of the classical dynamics 

The aim of this section is to explain how we numerically determine the classical dwell 
time Tn and the dimensionless scaling parameter rj appearing in Eq. fl58p that characterizes 
magnetic dephasing. To this end, we compute, with a ray-tracing algorithm, an ensemble 
of classical trajectories that enter the cavity from the left lead. The initial conditions 
of these trajectories are randomly selected from the intervals y^'^^ E [0, W] 

and p^^ G [— in a uniform manner, while we fix x^^^ = Xl and p° = \J p'^ — [p^y^Y (with 
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Figure D.18: Probability distributions (blue lines) for the path length L (left panel) and for the directed 
area A (right panel) for the case of the annular stadium billiard shown in Fig. m(b). Neglecting a short 
transient region, an exponential decay (red dashed lines) is fitted to these distributions. 



p = y/2mfi the total momentum of the classical particle). The propagation of a trajectory 
is continued until it exits the billiard via one of the leads. 

Fig. ID. 181 shows, for the case of the annular stadium billiard shown in Fig. [Il^(b), the 
numerically obtained probability distributions for the path length L and for the modulus of 
the directed area A that is accumulated along the trajectories according to Eq. (H^ . As it 



is expected for chaotic motion |46l. |83|. both probability distributions follow an exponential 
law after a short transient region. Fitting an exponential decay P{L) oc exp(— L/Lq) to the 
asymptotic behaviour of the probability distribution for the trajectory lengths, we obtain 
the dwell time via td = Lq/v with v ^ p/m the velocity of the particle. 

The distribution of directed areas P{A) can be determined from the distribution P{t, A) 
[see Eq. f l^ ] via 



P{A) 



P{t,A) exp{-t/TD)dt 



1 



: exp 



2\A\ 



(D.l) 



831 ]. The exponential decay 



and is also predicted to decrease exponentially with |^| |46l . 
P{A) oc exp (— l^l/^o) that is numerically encountered after a short transient region allows 
us to determine the characteristic scaling parameter via r] = 4:AI/{0,^^'^vt£i), using the dwell 
time Td that is obtained from the length distribution as explained above. 

Comparing the numerically computed dwell time rj^"™ -* with the "universal" prediction 



we obtain r 



(num.) 



D 



0.79{'kQ)/[{W + W)v]. This deviation is attributed to the finite 
width of the leads, which effectively compromises the assumption of ergodic motion that 
underlies the derivation of Eq. (H3|) . Regular islands of appreciable size, which might also 
give rise to a deviation of the dwell time from the universal prediction, could not be identified 
in the phase space of the two billiards. 
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Appendix E. Frequent integrals in the calculations of loop corrections 



Appendix E.l. The standard encounter integral 

In this appendix, we calculate the contribution of the encounter region. We first consider 
the absence of a nonlinearity. The corresponding integral is given by 

I^{E) = ( r ds r du——^-^^exp ( '-su] exp ( Jjd^] \ , (E.l) 



j:{E)tenc{su) \h J V ^' 

where l/Xfj=l/rD + n/rs accounts for the fact that we can have n G {0, 1, 2} stretches with 
a gauge field dependence within the encounter region, and = 27rmf2 is the volume of 

the energy shell in phase space. Following Ref. [60[, we first split the integration over u in 
two parts, 

/O /-c 
du . . . + du . . . , (E.2) 
-c io 

and make the variable transformation {s^u) {S,y) with S = sujc? G [—1,1] and 
y = c/\u\ = Tc/u G [1,1/5], with the associated Jacobian determinant c^/y, where the 
sign in the definition of y refers to the first and the second integral on the right-hand side of 
Eq. f lE.2p . respectively. In physical terms, we transform here from the phase space coordi- 
nates s, u to the action difference su measured in terms of c^, and to a coordinate y related 
to the time t„ = (1/A) ln(c/M) needed for the unstable phase space coordinate to evolve from 
the Poincare surface of section V to the limiting value =Fc. 

As tenci.su) = (1/A) ln(c^/|sM|) = tenc(|'S'|) docs uot dcpcud ou y, we obtain 

i/\s\ I / 1 X 

dy- = In { — ]= AWd^l) (E.3) 



1 y \\s 

for the integration over y. We then have 

-,2 \ /-l 



^"'^ dS cos ( ^] S'/^'^"^ 



As the limiting scale c for the coordinates s and u (i.e., the scale until which the linerization 
around the reference trajectory is still valid) generally depends on the energy the first term 
in Eq. flE.4p is expected to strongly oscillate when varying E and would thus vanish when 



performing the energy average. For the second term, we obtain after the transformation 

s' = Sc^/h 
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Figure E.19: Sketch of the two different possibihties for a nonhnearity event to enter the encounter region. In 
the scenario depicted on the left-hand side, the nonhnearity event moves along a stretch all the way through 
the encounter region. On the right-hand side, two of the four encounter stretches end at the nonhnearity 
block. This gives rise to a reduced encounter duration tenc{t,u) = t + tu which depends on the time interval 
t between the nonhnearity event and the Poincare surface of section V, and on the time t„ — (1/A) ln(c/|u|) 
between V and the borders of the encounter region. 



Assuming that the Ehrenfest time te = ln{c^/h) is much smaller than the dwell time 
td and the magnetic dephasing time scale tb, we have te <^ t„ as well as Ar^ ^ 1 and can 
approximate 

'y;^^l/(Ar„) 

exp 



Te 

Tr,. 



1 . 



(E.6) 



The remaining integral can then be evaluated in the semiclassical limit ^ — )■ by sending 
the upper limit of the integration flE.5|) to infinity, which finally yields 



I-{E) 



Ah 



dS 



sin {S) 
S 



2'nh 



n- 



(E.7) 



Appendix E.2. The encounter integral with an embedded nonhnearity event 

Now we consider the presence of a nonhnearity event within the encounter region. We 
first focus on the case that the nonhnearity event is moving along a stretch all the way 
through the encounter region, as depicted on the left-hand side of Fig. IE. 191 For this case, 
we have to evaluate the integral 




dsdu 



1 



S(E)tenc(sM) 



exp 



h^ 



su exp 



\su] 



(it exp 




where the additional integration variable t represents the location of the nonhnearity on 
a stretch within the encounter region. The gauge field dependence, manifested by the 
dephasing factor exp(— t/r^), emerges from the stretch along which the nonhnearity moves. 
We have 



(it exp 



t 

Tb 







^onc ( 


Tb 


1 — exp ^ 











(E.9) 
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which would also be obtained if the integrand in Eq. (1E.8|) was exp[— (tcnc — t) /tb], corre- 
sponding to the case that the other part of the stretch guiding the nonlinearity event would 
provide the gauge field dependence. Using the results from section Appendix E.l , we obtain 

1 



h{E) = rs [I'AE) - iliE)] 



(E.IO) 



We now analyze the second scenario, shown on the right-hand side of Fig. IE.19t where 
stretches of the encounter region terminate at a nonlinearity. The integral that has to be 
evaluated in this case is given by 




dsdu 



(l/A) ln(c/|s|) 



dt 



exp \ —su ] exp 




(E.ll) 

where we define tenc(t,u)=t + {1/ X)\n{c/\u\) as the reduced encounter time and l/r^ = 
1/td + u/tb, with n = 0, 1, 2 the number of pairs of imbalanced stretches that give rise to a 
gauge field dependence. As indicated in Fig. IE.19t the integration variable t represents the 
time between the nonlinearity and the Poincare surface of section V within which the stable 
and unstable coordinates are defined. 

Following Refs. 60|, |6l|, |62|, we split, as in Section Appendix E.l, the integration over 
u according to Eq. f lE.2p and make the variable transformation (s, u, t) i— )■ {S, y, i) with 
S = su/c^ e [-1,1], t = i^nc{t,u) =t+(l/A)ln(c/|M|) e [0, (l/A)ln(l/|5|)], and y = c/\u\ G 
[l,exp(Af)], with the associated Jacobian determinant c^/y. The integration over y yields 



/•exp(Af) 1 

/ dy- 

Ji y 



xt. 



(E.12) 



We then evaluate 



2c^X 
'2c^Ar„ 
'AfXTn 



dS 



(l/A) ln{l/|SI) 



(it exp 

dS [1 - ISl^/^^-^^ exp 



-1 

h / c 

— sm — 



dS cos 




(E.13) 



The first contribution in the last line of Eq. (IE.13P vanishes when performing the energy 
average, whereas the second term yields [—{hn)/{2Xc^Tn)], as seen in Section Appendix E.l 
We thus obtain 
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